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ABSTRACT 

Rossby waves (r-modes) in rapidly rotating neutron stars are unstable because of the emis- 
sion of gravitational radiation. As a result, the stellar rotational energy is converted into both 
gravitational waves and r-mode energy. The saturation level for the r-mode energy is a funda- 
mental parameter needed to determine how fast the neutron star spins down, as well as whether 
gravitational waves will be detectable. In this paper, we study saturation by nonlinear trans- 
fer of energy to the sea of stellar "inertial" oscillation modes which arise in rotating stars with 
negligible buoyancy and elastic restoring forces. 

We present detailed calculations of stellar inertial modes in the WKB limit, their linear 
damping by bulk and shear viscosity, and the nonlinear coupling forces among these modes. The 
saturation amplitude is derived in the extreme limits of strong or weak driving by radiation 
reaction, as compared to the damping rate of low order inertial modes. In the weak driving case, 
energy can be stably transferred to a small number of modes, which damp the energy as heat or 
neutrinos. In the strong driving case, we show that a turbulent cascade develops, with a constant 
flux of energy to large wavenumbers and small frequencies where it is damped by shear viscosity. 

We find the saturation energy is extremely small, at least four orders of magnitude smaller 
than that found by previous investigators. We show that the large saturation energy found in 
the simulations of Lindblom et al. (2001, 2002) is an artifact of their unphysically large radiation 
reaction force. In most physical situations of interest, for either nascent, rapidly rotating neutron 
stars, or neutron stars being spun up by accretion in Low Mass X-ray Binaries (LMXB's), the 
strong driving limit is appropriate and the saturation energy is roughly £' r _ mo d e /(0.5Af r^fl 2 ) ~ 
0A-fg r /il ~ 10~ 6 (i4pi n /10 3 Hz) 5 , where M and r» are the stellar mass and radius, j gr is the 
driving rate by gravitational radiation, f2 is the angular velocity of the star, and i4 P i n is the 
spin frequency. At such a low saturation amplitude, the characteristic time for the star to exit 
the region of r-mode instability is >10 3-4 years, depending sensitively on the instability curve. 
Although our saturation amplitude is smaller than that found by previous investigators, it is still 
sufficiently large to explain the observed period clustering in LMXB's. We find that the r-mode 
signal from both newly born neutron stars and LMXB's in the spin down phase of Levin's limit 
cycle will be detectable by enhanced LIGO detectors out to ~ 100 — 200 kpc. 

Subject headings: stars: neutron — gravitational waves — turbulence — stars: oscillations 

Introduction Theoretically, we expect neutron stars can ro- 

tate up to ~ 10 3 Hz without breaking apart (Cook 
What sets the observed spin rates of neutron et aL 1994^. Fryer & H eger 2000; Heger et al. 

11 s ^ 2000). However, for the rapidly accreting, weakly 



magnetic LMXB's, oscillations seen during type 
I X-ray bursts (Strohmayer et al. 1996), as well 
as quasi-periodic oscillations (van der Klis 1998), 
seem to indicate spin frequencies narrowly clus- 
tered near 300 Hz. If LMXB's are the progeni- 
tors of millisecond pulsars, and the timcscale over 
which they should be spun up by accretion is only 
<~ 10 7 yr for high accretion rates, why aren't more 
stars spun up near 10 3 Hz over their > 10 9 year 
lifetime? 

Wagoner (1984) proposed that for weakly mag- 
netic neutron stars, the spin up torque due to 
accretion is balanced by spin down torque from 
gravitational radiation reaction. There are cur- 
rently two distinct models to explain the non- 
axisymmetric deformation of the star producing 
the radiation. The first mechanism involves mass 
quadrupole deformations of the neutron star crust 
(Bildsten 1998; Ushomirsky et al. 2000) while the 
second involves mass-current quadrupole emission 
from the r-mode instability (Bildsten 1998; Ander- 
sson et al. 1999b, 2000), which will be examined 
in detail in this paper. 

Many young neutron stars associated with su- 
pernova remnants also seem to be spinning slowly, 
in spite of the theoretical expectation (Fryer & 
Heger 2000; Heger et al. 2000) that typical 8 - 
25 Mq progenitors lead to neutron stars rotat- 
ing with periods of order 1 ~ 1 msec. Kaspi 
& Helfand (2002) cite the following examples for 
the inferred initial spin period P n i t and age T of 
some of the fastest rotators: the Crab pulsar with 
P init = 19 msec and T = 948 yr; PSR J0537- 
6910 in host remnant N157B with P; n i t < 10 msec 
and T = 5000 yr; PSR B1951+32 in CTB 80 has 
Pnit < 39 msec and T = 10 5 yr. The Crab is 
by far the most certain estimate for P; n it, with a 
known age from the historical supernova and mea- 
sured braking index. However, Kaspi & Helfand 
(2002) also note several slow rotators, such as PSR 
J1811-1925 in Gil. 2-0.3 with P init = 62 msec and 
age T = 2000 yr. 

The apparent discrepancy between the theoreti- 
cally expected fast rotation rates and the observed 
slow rotation rates could be reconciled if some 



1 This result depends sensitively on the poorly understood 
coupling between the core and envelope of the progenitor. 
Angular momentum transport mechanisms due to, for in- 
stance, weak magnetic fields may decrease the rotation rate 
of the core prior to collapse. 



mechanism could slow down fast rotators, effec- 
tively preventing them from reaching millisecond 
spin rates. The r-mode instability is a possible 
mechanism. 

This instability was discovered by Andersson 
(1998) and Friedman & Morsink (1998) showed 
that all rotating, inviscid stars are unstable be- 
cause of this general relativistic effect. The insta- 
bility arises when certain stellar oscillation modes, 
called Rossby waves (or r- modes), are driven un- 
stable by the emission of gravitational waves. As 
a result, the rotational energy of the star is con- 
verted into both mode energy and gravitational 
waves, causing the star to spin down. Detailed cal- 
culations (Lindblom, Owen & Morsink 1998; An- 
dersson et al. 1999a; Kokkotas & Stergioulas 1999; 
Lindblom et al. 1999; Bildsten & Ushomirsky 
2000; Levin & Ushomirsky 2001; Lockitch & Fried- 
man 1999) show that viscous dissipation by large 
scale shear, boundary layer shear at the crust-core 
interface, and modified URCA bulk viscosity are 
likely insufficient to counter this driving in rapidly 
rotating neutron stars. However, Lindblom & 
Owen (2002a) point out an interesting mechanism 
for bulk viscosity arising from hyperon interac- 
tions which may overcome the driving. Mcndell 
(2001) has investigated the effects of magnetic 
fields on the boundary layer, finding that large 
fields can significantly increase the damping rate. 
Lastly, the work of Levin & Ushomirsky (2001) 
shows that damping from the crust-core boundary 
layer leads to a double-valued instability curve, 
which may explain why LMXB spin frequencies 
are lower than those of the millisecond pulsars. 

The instability may be important in two re- 
spects. First, r-modes in any neutron star ro- 
tating faster than some critical rate will become 
unstable, causing the star to rapidly spin down. 
Hence, r-modes may set a maximum rotation rate 
for neutron stars. Second, the enormous amount 
of energy radiated in gravitational waves may be 
detectable by LIGO. 

In section 2 we review how nonlinear satura- 
tion occurs in the limits of weak and strong driv- 
ing. We derive formal expressions for the satura- 
tion amplitude, which depend on the microphysi- 
cal details of the nonlinear interaction and damp- 
ing rates. Section 2.1 contains a review of nonlin- 
ear coupling of just three oscillation modes, with 
emphasis on amplitude saturation by the paramet- 
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ric instability. Section 2.2 reviews amplitude sat- 
uration by a continuum of modes in which a well 
defined inertial range exists. In section 3, we dis- 
cuss the modes present in rapidly rotating neu- 
tron stars, arguing that the buoyancy and elastic 
restoring forces are weak compared to the Coriolis 
force. We compute WKB inertial cigenmodes in 
section 3.2. The nonlinear coupling coefficients are 
computed in section 4, and damping rates in sec- 
tion 5. The saturation amplitude for the discrete 
case is discussed in section 6, and the continuum 
case in section 7. Neutron star spin evolution due 
to the r-modc instability is discussed in section 9. 
Our results are compared to those of previous in- 
vestigators in section 8. We discuss the detectabil- 
ity of the gravitational wave signal in section 10, 
and give a summary and conclusions in section 1 1 . 
Two appendices give detailed calculations of the 
turbulent cascade for stellar inertial modes, and 
the nonlinear force coefficients. 

2. Saturation by Nonlinear Mode Cou- 
pling 

We start by reviewing the equations of motion 
for the mode amplitudes, and then specialize to 
the weak and strong driving limits. 

We will work in a reference frame co-rotating 
with the star. Expansion of the fluid displace- 
ments, relative to the co-rotating frame, in terms 
of the linear eigenmodes 
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leads to the following system of coupled oscilla- 
tor equations for the dimcnsionless complex am- 
plitudes q a (t) (Schenk et al. 2002): 



q a + iuj a q a 
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Qfa (2) 



The left hand side of cq.2 represents an unforced 
oscillator of rotating frame frequency oj a , while 
the terms on the right hand side are the driving 
(+) or damping (— ) term and the nonlinear term, 
which is quadratic in q. In our notation, K a( g 7 is 
roughly the ratio of interaction energy to mode en- 
ergy at unit amplitude. The rotating frame mode 
energy is E = 2£ , un j t |g| 2 , where £7 U nit is a (arbi- 
trary) unit of energy which we find convenient to 



set to -E un it = 0.5Mr 2 fl 2 2 . Here M and r* are the 
stellar mass and radius, and fi = fle z is the an- 
gular velocity. The sum over modes involves 
a sum over the mode {ojp, $,p), with amplitude qp, 
as well as its complex conjugate (— ujp,£*p), with 
amplitude q% (see Schenk et al. 2002 for a detailed 
derivation, but note that the type of index denoted 
there by A is denoted here by a.) 

2.1. the discrete limit 

In the regime where the driving rate of the un- 
stable mode is smaller than the damping rates of 
low order modes, the instability can be saturated 
by a transferal of energy to a small number of 
damped modes. We will begin by discussing the 
coupling between the "parent" r-mode, and two 
damped "daughter" modes. Although an ideal- 
ization, this basic problem is soluble, and indi- 
cates which modes couple most strongly to the r- 
mode. We review the dynamics of such 3-modc 
networks, including the parametric instability, the 
equilibrium solution, and the linear and nonlin- 
ear stability of the equilibrium solution (Wersinger 
et al. 1980; Wu & Goldreich 2001; Dziembowski & 
Krolikowska 1985; Dimant 2000; Abarbanel et al. 
1993). 

In terms of the real amplitude and phase vari- 
ables, defined by q — Aexp(-itp), the equations 
for a system of three modes are 



Ai = 7i-Ai — ujikA 2 As simp 
A 2 = -72^2 - lu 2 kA 3 Ai sin tp 
A3 = -73^3 - uznA x A 2 sin tp 



tp = 5lu — k cos ipAiA 2 A 3 



ui 2 
AV A\ 



A 3 



where the index 1 refers to the parent and 2 and 3 
refer to the two daughter modes. We have defined 
the frequency detuning 5u> = u)\ + lo 2 + 0^3 , cou- 
pling coefficients K123 = «exp(— id), and the rela- 
tive phase tp = 8 + pi + <p 2 + tp 3 . The qualitative 
features of the time evolution, such as equilibrium 
and stability, depend only on the three dimension- 
less parameters (72 +73V71, 72/73, and 5w/7i- 

The parametric instability (Landau & Lifshitz 
1969; Dziembowski & Krolikowska 1985; Kumar 
& Goodman 1996; Wu & Goldreich 2001) is a 



2 The nonlinear interaction energy also scales as the rotation 
energy of the star 
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Fig. 1. — Saturation in the discrete limit. Mode amplitudes are evolved in time for the case where 72/71 = 
73/71 = 3, and 800/^1 — 10. The solid line is the driven mode, and the dashed line represents the two 
daughter modes, which are identical in this example. The absolute scale is arbitrary. The amplitudes and 
phase are started well below their equilibrium values, which are denoted by the straight upper solid line and 
the straight dashed line. The lower straight solid line is the parametric threshold. Notice that the daughter 
mode amplitude decreases until the parent exceeds the threshold, at which point the daughter amplitude 
rises exponentially. Qualitatively similar evolutions are obtained for a wide variety of initial conditions. The 
parameters 71, 72, 73 and Su> were chosen so that the solution would converge to the equilibrium values. 



mechanism by which the daughter mode ampli- 
tudes will grow exponentially if the parent mode 
amplitude exceeds a certain threshold 3 . The re- 
sult is that energy can be quickly taken out of 
the parent mode and transferred to the daugh- 
ter modes, providing a means to limit the parent 
mode's amplitude. Furthermore, the growth rate 
of the daughters is larger than the growth rate of 
the parent so that the daughters can "catch up" 
to the parent even if they start from a lower am- 
plitude. 

The parametric instability can be derived in the 
approximation where the parent mode's amplitude 
is much larger than the daughter modes' ampli- 
tudes so that the influence of the daughters on the 
parent can be neglected. Performing a linear sta- 
bility analysis of eq.3 (Landau & Lifshitz 1969), 
one finds that the daughters grow exponentially 
like exp(i/r) when the parent mode amplitude A\ 



3 A simple example of parametric instability is a pendulum 
in which the length of the string is being varied periodically. 
See Landau & Lifshitz (1969). 



exceeds the critical value given by 



Ai = 



7273 

K 2 Ul2U>3 
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72 + 73 



(4) 



where 72,3 = 72,3 + 1/r. In particular, the thresh- 
old at which the instability first starts to operate 
is given by eq. (4) at r = 00, 



A\ = 



k 2 Q 2 Q3 



1 + 



5u> 



72 +73 



(5) 



where Q2 — ^2/72 and Q3 = W3/73 are the qual- 
ity factors of the daughter modes. In the limit 
72 j 73 - * of negligible damping of the daugh- 
ter modes, it is useful to consider in addition 
the threshold above which the daughter mode's 
growth rate will exceed that of the parent mode. 
This is given by eq. (4) at t = 7f , 72 = 73 = 0: 

We give an example showing the parametric insta- 
bility in Fig. 1. 
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Once the parametric instability occurs, the 
daughter modes start to grow rapidly. We now 
discuss the conditions under which the subsequent 
evolution leads to a saturation of the parent mode 
in the three mode system. 

Setting the time derivatives in eq.3 to zero, one 
finds the equilibrium solution for the parent 



A 2 = 
1 k 2 Q 2 Q 3 



1 
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72 + 73 - 7i 
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and daughter mode energies A\^jA\ = Qi,zlQ\-, 
where Q — uj/j is the quality factor of the mode. 
Naively, one expects that energy transfer from the 
parent to the daughters occurs only if the daugh- 
ter modes have a lower energy than the parent, 
implying a lower daughter mode quality factor. 
This expectation is verified by a stability analy- 
sis (Wu & Goldrcich 2001; Dimant 2000) which 
shows that the equilibrium solution is stable only 
when 72 + 73 > 71 is satisfied. 

More precisely, there are three different regimes 
in the three dimensional space of parameters (72 + 
73) /7ii 72/73 and 5u>/ji. First, the equilibrium 
solution is linearly stable to small perturbations 
if two conditions are met (Wersinger et al. 1980): 
(i) the ratio of damping to driving is sufficiently 
large 72 + 73 ^ 71, and (ii) the detuning is suf- 
ficiently large, Slo > (72 + 73)/2. Second, in the 
regime where 72 + 73 ^ 71 but where the detun- 
ing is small, Suj < (72 + 73)/2, the the amplitudes 
and phase undergo limit cycles characterized by 
bounded, quasiperiodic orbits, as shown by Fig. 2 
of Wu & Goldreich (2001). Those limit cycle solu- 
tions have time averaged parent mode amplitudes 
comparable to the equilibrium amplitude (7), so 
the equilibrium amplitude still characterizes the 
motion. Third, if the daughter mode damping is 
insufficient, 72 + 73 < 71, all three amplitudes rise 
without bound and the solution is nonlinearly un- 
stable (Dimant 2000). For our purposes, any so- 
lution which is nonlinearly stable can saturate the 
growth of the r-mode, so that the effective stability 
criterion is 



72 + 73 ^71- 



(8) 



In the regime where the equilibrium solution 
is stable, it acts like an attractor, and the sys- 
tem tends to evolve into this equilibrium after the 
daughter mode amplitudes become comparable to 



that of the parent mode. The example shown in 
Fig. 1 exhibits this behavior, even though the sys- 
tem is started well away from equilibrium. Note 
that the equilibrium parent mode amplitude (7) is 
always approximately equal to the threshold am- 
plitude (5), in the regime (8) where the energy 
transfer is stable. 

The parametric instability can provide a means 
for saturating the r-mode amplitude. Suppose 
that a daughter pair exists which is parametri- 
cally unstable for a certain value A\ of the par- 
ent mode amplitude, and that no other daughter 
pairs are unstable at that amplitude. Then, if the 
resonance is sharp, it is plausible that only the 
parent and two daughter modes are relevant, and 
if the condition (8) is satisfied so that the transfer 
of energy is stable, then driving of the r-mode by 
gravitational radiation reaction can be balanced 
by nonlinear energy transfer to the pair of daugh- 
ter modes. Thus, the daughter mode pair for which 
the instability threshold (5) is lowest sets the sat- 
uration amplitude for the r-mode, if the stability 
constraint (8) is satisfied for that daughter mode 
pair. Daughter pairs with higher thresholds will 
not be excited because the parent's amplitude can- 
not rise much above the lowest threshold (see Fig. 
1)- 

The task of finding the saturation amplitude 
in the weak driving regime involves searching 
through all possible daughter mode pairs to mini- 
mize the parametric threshold (5). This amounts 
to maximizing n while minimizing the mismatch 
Suj 2 + (72 +73) 2 subject to the stability constraint. 
Once this "best" daughter mode pair has been 
found, the saturation amplitude is 
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I best pair, 



(9) 



assuming that the strong resonance condition Su>< 
72+73 is satisfied. We can state the following rule 
of thumb: for coupling coefficients of order unity, 
the r-mode will saturate to an amplitude less than 
unity if the best daughter pair are high Q oscilla- 
tors. Quality factors of low order global modes in 
neutron stars can easily be 10 6 or larger. 

Finding the saturation amplitude in the weak 
driving regime has now been reduced to the fol- 
lowing physics problem. First determine the oscil- 
lation modes present in the star. Calculate their 
damping and driving rates, as well as the nonlin- 
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car coupling coefficients between daughter pairs 
and the r-mode. Once the magnitude and scalings 
of these quantities are known, reliable estimates of 
the parametric threshold can be made (see Sec. 6 
below) . 

Finally, note that nonlinear coupling terms such 
as (72 — K2ii5i which couple the parent mode 
twice with a daughter mode have been ignored in 
cq.3. Since these terms scale as A\, instead of 
A\ as for the parametrically driven modes, they 
are smaller in the weakly nonlinear regime. In ad- 
dition, the coupling coefficients drop off rapidly 
for this type of coupling as the wavenumber of 
mode 2 is increased (see appendix B). Hence, 
only daughter modes with comparable wavenum- 
ber to the parent couple well. However, the res- 
onance condition cannot be finely tuned for com- 
parable wavenumber modes, since there are so few 
of them. As opposed to the couplings K2ii<Zi, 
parametric type couplings have the double advan- 
tage of allowing coupling of the parent mode with 
daughter modes of arbitrarily large wavenumber, 
and the resonance condition becomes satisfied to a 
higher degree of precision for large daughter mode 
wavenumber. 

2.2. the continuum limit 

In the above "discrete" scenario, the saturation 
amplitude of the driven mode scales as A ~ Q^ 1 , 
where Qd is the quality factor of a damped mode. 
In the "continuum" picture that we now discuss, 
the saturation amplitude is independent of the lin- 
ear damping rates, since the energy is transferred 
by nonlinear interactions. In this cascade picture, 
both the shape and normalization of the wave en- 
ergy spectrum are set only by the detailed non- 
linear interaction between waves, and the energy 
input to the system. 

How does the cascade arise? Imagine starting 
with a system in the weak driving limit and adi- 
abatically increasing the driving. When j gr be- 
comes greater than the damping rate 72 + 73 of 
the daughter pair with the lowest threshold, the 
equilibrium solution for that pair is no longer sta- 
ble and the energy of all three modes will begin 
to grow. When the energy has grown to the point 
that additional parametric thresholds are crossed, 
and if the energy transfer to these pairs is stable, 
the driven mode will again be saturated. As the 
driving is increased, this process will continue until 



many daughter modes are excited with large am- 
plitude. Since linear damping is smaller than driv- 
ing over a certain range of daughter mode length- 
scales, an inertia! range has formed where nonlin- 
ear forces are dominant. We now proceed to give 
a heuristic derivation of instability saturation in 
this continuum case, leaving the detailed deriva- 
tion appendix A. 

Since many modes are excited, we treat the 
quantum numbers for each mode as a continuum. 
Introducing the "occupation number" 4 (quasi- 
particle number) for mode a 

JVa = |«a|VM > 0, (10) 

the mode energy becomes 

E a — 2_E unit |<? Q | 2 = 2E un it\uJ a \N a . (11) 

Eq.2 describes both the fast variation of each in- 
dividual mode, as well as the slow variation due 
to nonlinear interactions between modes. We may 
average over the fast oscillations using the ran- 
dom phase approximation (Zakharov et al. 1992; 
Kumar & Goldrcich 1989; Wu 1998) if the phase 
randomization time set by the wave dispersion is 
shorter than the nonlinear interaction timescale. 
Since the dispersion time is comparable to the 
mode period for inertial waves, this is equivalent 
to the weak nonlinearity condition. The resultant 
kinetic equation for the wave amplitudes is (Za- 
kharov et al. 1992; Kumar & Goldreich 1989) 

N a = I a +T a N a (12) 

where I a represents the rate of change of N a due 
to nonlinear interactions, and T a is the rate of 
driving (> 0) or damping (< 0). The interaction 
term has the form 

(3-f 

x (saNpN^ + s p N 7 N a + s 7 N a N p ) (13) 

where s a is the sign of the frequency of mode a. 

To proceed further, we must introduce a few 
properties of the oscillation modes to be derived 
in section 3. Let n, k, and m denote the perpen- 
dicular (to ft), parallel, and azimuthal number of 

4 To get a quantity with the units of action, multiply N a by 

^unit- 
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Fig. 2. — Saturation in the continuum limit. Energy is input to the r-mode at the outer scale of both 
wavenumber n and frequency /i. Energy then cascades to small \x and large n along the integral curves 
n a„^a n _ (ooust) of the energy flux vector field (16). 
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nodes, respectively. Since inertial mode oscillation 
frequencies are proportional to the rotation fre- 
quency we write the mode frequency asw = 2£l^i, 
where /j. 2 < 1 is the dimensionless frequency. 

Approximate stationary solutions of eq.12 are 
found in two steps (see Zakharov et al. (1992) for 
detailed derivations). First, one ignores the driv- 
ing and damping, so that the energy flux is con- 
served. In this case, the energy flux T is defined 

by 

w J dml = -V k -F (14) 

where is the gradient in momentum space and 
we have integrated over the m quantum number. 
In appendix A, we show that stellar inertial waves 
support a flux of energy to large n and small \i. A 
schematic drawing of this cascade is given in fig. 2. 



The occupation number 5 for each mode is 

N = n- 1 JV n-'Vr 1/2 oc n-y 2 \k\-^ 2 {15) 

where the normalization constant is related to the 
fluxes T n and in the n and directions, re- 
spectively, by 

T n = S^^n-Vr'^unit 

T» = -8a^n- 2 ^nE unit . (16) 

The constants a n and are order unity and pos- 
itive. A surface of constant energy in momen- 
tum space has fi cx n 8 , showing that the energy 
cascades to small frequencies quite rapidly with 
wavenumber, because of the strong dependence 
of the coupling coefficients on frequency (see ap- 
pendix B). 

The final step is to match the inertial range so- 
lution to the driving range. In other words, we 

5 The scaling of this expression for the occupation number 
can be simply derived from dimensional analysis together 
with the fact that 3-mode interactions dominate over 4- 
modc and higher order interactions. See, e.g., Sec. 3.3.1 of 
Zakharov et al. (1992). 
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need to find steady state solutions to the equation 
I + IgrN = 0, where -f gr is the driving rate by 
gravitational radiation. We approximate N in the 
driving region by extending the inertial range so- 
lution. The power input by the instability is given 
by 6 

input power = J dndk dm^ gr E 
2 8 

~ -n 3 r -/ gr E r = -TgrSunitiVoVVr 72 (17) 

where the r subscript denotes the driven r-modc 
and we have approximated the (dimensionless) 
volume in phase space being driven as oc nf. The 
energy escaping the driving region is given by 
integrating up the flux through each boundary. 
Roughly, this is given by 



/■Mr 

output power = dfi nT n 

•V(l-l/n r ) 
rn r +l 

+ dn n\^\ 

J n r 

,-1 Ar2/ 



Equating the input to output power we can solve 
for the normalization constant 



1 



1 



3a n + ft 



1ml,, 1 / 2 

H'r ' 



(19) 



Given the normalization, we find 



E n 



E(r — mode) 
0.5Mr 2 Q 2 



3(c 



a, 



(20) 



where we parametrize our inexact treatment of 
the matching condition with the parameter a e — 
2n~ 4 (i r /3(a n + a M ). If we use the quantum num- 
bers of the r-mode, n r = 3 and [i r = 1/3, we 
find a e ~ 4 x 10~ 4 . However, we choose to be 
very cautious about this factor since we are ex- 
trapolating a WKB treatment into the regime of 



6 Even if a relatively narrow region ~ An 3 in phase space is 
being driven, Zakharov et al. (1992) find that one should 
use the whole volume ~ instead of An 3 since a peak 
develops in the driving region. Since the r-mode has rel- 
atively small wavenumber, the width of the driving region 
occupied by the r-mode may be considered relatively wide 
(An/n ~ 1/3, 8k/k ~ 1, Sm/m ~ 1/2). 



low order modes 7 . A more conservative estimate 
would be to set n r = fi r = 1, giving a e — 0.1. We 
will use the more conservative result for numerical 
work in the rest of this paper, but recall that it 
may overestimate the saturation amplitude by up 
to three orders of magnitude. Using the r-mode 
driving rate from eq. 58, the final result for the 
saturation energy is 



E{r — mode) 
Q.hMrlVL 2 



- (^) (21) 



where ^khz is the spin frequency in units of 
1000 Hz. 

Why is the saturation amplitude so small? The 
factor 7 9r /fi is inevitable 8 since the only quantity 
with the units of frequency in the nonlinear inter- 
action rate is fL The numerical factor a e depends 
on considerations such as the effective volume and 
area of the driving region, and the power in the 
driving region relative to the largest scale (energy 
bearing) waves. 

Eq.21 is one of the central results of this paper. 
It applies when nonlinear energy transfer is faster 
than linear damping. If nonlinear energy transfer 
becomes slower than linear damping, the discrete 
limit of section 2.1 is recovered. Note that the 
saturation amplitude decreases very rapidly with 
stellar spin frequency. 

3. Oscillation modes in rapidly rotating 
neutron stars 

In this section we discuss the oscillation modes 
present in rapidly rotating neutron stars. We ar- 
gue that at the rapid rotation rates of interest for 
the r-mode instability, the buoyancy and elastic 
restoring forces can be ignored in comparison with 
the Coriolis force. The resulting modes which are 
restored by the Coriolis force are called inertial 
modes, of which the r-modes are a subset. 



7 The detuning may become non-negligible for low order 
modes. In addition, the resonance width from 7 gr can 
become important for coupling directly to the unstable r- 
mode. 

8 If the energy transfer is local in frequency space, this scaling 
will also hold for interaction with other wave families, such 
as inertial-gravity modes. 
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3.1. motivation for inertial waves 

Within a minute after their birth in a super- 
nova, neutron star cores have become transpar- 
ent to neutrinos and cooled down sufficiently to 
form a degenerate gas of neutrons, with a small 
admixture of electrons and protons determined by 
beta equilibrium. As shown clearly by Reiseneg- 
ger & Goldreich (1992), the varying electron frac- 
tion y e ~ 6.0 x lCT 3 0/2.8 x 10 14 g cm" 3 ) in 
the star causes a stable stratification and result- 
ing buoyancy force: Since the neutron pressure 
p n oc [p/(l + j/ e )] 5 / 3 , displacing a fluid element up- 
ward on timescales slower than the sound crossing 
time and faster than the timescale of the beta re- 
actions results in the fluid element being heavier 
than its surroundings, since it came from a region 
of larger y e . An oscillatory motion results, with 
maximum frequency of order the Brunt- Vaisala 
frequency (Reisenegger & Goldreich 1992) Nb v ~ 
{O.Syeg/Hp) 1 / 2 ~ 500 sec" 1 - 2?r x 100 Hz, where 
g is the local gravity and H p is the local pressure 
scale height. 

The buoyancy force on a fluid element is just 
-Fbuoy = —gdp — —pN£ v l; r , where £ r is the radial 
component of the Lagrangian fluid displacement. 
The Coriolis force is given by F cor <~ 2pQw£, so 
that the ratio of these two forces for u <~ fi is 
roughly 




In a detailed study of the solutions to the fluid 
perturbation equations for rotating stars includ- 
ing buoyancy, Yoshida & Lee (2000) showed that 
in the limit of Fb uoy /F COI » 1 the solutions are 
approximated very well by the r- and g-modes. 
This limit of large buoyancy force was examined 
by Morsink (2002) who showed that the nonlinear 
couplings between r-modes are too small to cause 
saturation to occur. In the limit of -Fbuoy/-Pcor <C 1 
Yoshida & Lee (2000) have shown that the solu- 
tions of the perturbation equations are well ap- 
proximated by the inertial modes. As long as we 
restrict our calculations to stars spinning at a fre- 
quency greater than 100 Hz, the inertial modes 
with zero buoyancy are very good approximation. 
As we are interested in the possibility of mode 
saturation at spin frequencies at least as large as 
300 Hz, the inertial modes are the most relevant 



modes and it is possible for us to ignore the buoy- 
ancy force as a first approximation. This enables 
us to find simple solutions for the modes if we fur- 
ther approximate the shape of the star as spheri- 
cal, a valid approximation for rotation rates well 
below the breakup rate. However, we expect that 
the qualitative results found here will hold even 
in the case when buoyancy is included. The rea- 
son is that the approximations made still provide 
a dense spectrum of modes that may be arbitrarily 
resonant with the r-mode in the continuum limit. 
Although the numerical value of the coupling coef- 
ficients and damping rates may change because we 
don't have exactly the correct shape of the eigen- 
functions, we are confident that the essential qual- 
itative features present in our simple example will 
carry over. 

Levin & Ushomirsky (2001) have shown that 
the clastic restoring force in the neutron star 
crust becomes small compared to the Coriolis force 
above a rotation rate of <~ 50 Hz. The net result is 
that core modes can penetrate into the crust, with 
only a small discontinuity at the crust-core bound- 
ary because of the impedance mismatch. We will 
ignore crustal elasticity for the remainder of this 
paper. 

We have not included superfluid effects in our 
calculations. The principal new effect would be 
dissipation due to mutual friction (the modes 
themselves are not expected to be changed very 
much; see, e.g., Lindblom & Mendell (2000)). 
However, we note that our estimate of the sat- 
uration amplitude does not depend on the dissi- 
pation rate if an inertial cascade forms (although 
the outer scale of the inertial range might be af- 
fected). Thus, if saturation involves a cascade of 
energy to numerous inertial modes, we still ex- 
pect our estimates to hold. Our estimates would 
change if dissipation via mutual friction is strong 
enough that only a few modes are excited para- 
metrically. We postpone a thorough examination 
of this case, which would depend on uncertain mu- 
tual friction coefficients, for another paper. How- 
ever, either way, the saturation amplitude will still 
be very small. 

In the next subsection, we discuss inertial mode 
eigenfunctions in weakly stratified stars. 
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3.2. stellar inertial modes 

We solve the Euler and continuity equations for 
adiabatic perturbations of a rotating star. The 
background star is taken to be spherically symmet- 
ric with uniform rotation rate = fle z and neg- 
ligible stable stratification. Perturbation modes 
of the form cxp(im<f) — iut) are found using the 
Cowling approximation. 

The Euler, continuity, and state equations are 



pi + 2pfl x £ 



-VSp + gSp (23) 



Sp = c 2 Sp, 



(24) 



(25) 



where we have ignored terms involving the Brunt - 
Vaisala frequency 



N bv = 



dr 



1/2 



(26) 



a valid assumption for u> ^ Nb v and fl ^ Nf, v . 
Here we have defined the adiabatic index F\. In 
this limit, the adiabatic sound speed c and density 
scale height H are related by c ~ {gH) 1 / 2 . Sub- 
stituting the assumed dependence on cf) and t, and 
eliminating Sp, we find that eqns. 23 - 25 become 



£ + iq e z x £ = 

LU 2 C 

c l H 



(27) 
(28) 



Here we have replaced the Eulcrian pressure per- 
turbation by the quantity ip defined by Sp = pco 2 tpi 
and defined the dimcnsionless inverse frequency 
q = 20/ 'u>. We will also heavily use the dimen- 
sionless frequency p = l/q = uj/20,. We drop the 
term aj 2, ip/c 2 oc fl 2 /(GM/rl) since we are work- 
ing to leading order in Q; this is consistent with 
our assumption that the background star is spher- 
ical and suffices to compute the mode functions to 
leading order in f2. 

The Euler equation 27 can be solved 9 for £ in 
terms of -0: 

£ = (l-q 2 )- 1 [ViP-q 2 e z (e z -Vi;)-iqe z : 



' The determinant of this transformation is singular only if 
io 2 =4H 2 . 



Substituting eq. 29 into the continuity equation 
28 gives the wave equation 



^2 , 2^-0 



q cos( 

or 



dz 



mq 



The boundary condition near the surface is 
that the Lagrangian change in the pressure van- 
ish, Ap = Sp + £ ■ Vp = 0, so that Sp ~ pgC ■ 
This is just the statement that the surface layer is 
hydrostatic, a consequence of the vanishing sound 
crossing time across a scale height for small depth. 

Equation 30 does not appear to be solvable by 
separation of variables 10 . This motivates us to ex- 
amine approximate solutions valid for short wave- 
lengths. Our solution generalizes the exact solu- 
tion of Bryan (1889) for the constant density star; 
in fact our solution is just Bryan's solution divided 

by VP- 
Defining 



ip(x) 



1/2 

^) /(*) 
P 



(31) 



where ipo is a normalization constant and po is the 
central density, the differential equation for / is 



V 2 / 



dz 2 



JC 2 f 







(32) 



with 

K 2 



(l-gW0) - 



ldH~ 



+ 



mq 



2rH 



2 dr 



q 2 sin 2 i 



4H 2 



(33) 



The first two terms in eq. 32 are just the usual dif- 
ferential equation for inertial modes, as derived by 
Bryan. The compression term £, r /H in the conti- 
nuity equation 28 is imaginary in the WKB limit, 
and leads to the WKB envelope p -1 / 2 . The defini- 
tion in eq.31 accounts for this envelope, so that the 
correction terms in eq.33 are now real. The K? f 
term is most important near the surface, where it 
scales as (\k\H)~ 2 relative to the other terms in 
cq.32 (|fe| is the local WKB wavenumber). This 



10 By separation of variables, we mean that (1) the differen- 
tial equation is separable, and (2) the boundary conditions 
are applied on a surface where one of the coordinates is 
constant. 
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term describes a slow variation of the wavenum- 
ber with position, and is negligible for short wave- 
length modes. Henceforth, we set JC ~ in the 
interior of the star. 

The short wavelength approximation breaks 
down when the vertical wavenumber k r 11 is com- 
parable to the scale height, at about 2k r H ~ 1, 
as can be verified directly from eq. 32. When 
the eigenfunction is constant over a scale height 
one may picture that the wave attempts to lift an 
entire scale height of material, which causes re- 
flection. We extend our interior solution to the 
surface by replacing p in eq. 31 with a "cut- 
off" value p cut which becomes constant about one 
wavelength from the surface. 12 

Bryan (1889) found a solution to cq.32 for 
JC = in terms of an ingenious bi-sphcroidal co- 
ordinate system that depends on the frequency p. 
Lindblom & Ipser (1999) have given a careful dis- 
cussion of this coordinate system, paying partic- 
ular attention to the behavior of the coordinates 
at the surface of the star. Define the bi-spheroidal 
coordinates \p\ < x\ < 1 and — \p\ < x 2 < \p\ by 



r* 



y = 



z = 



1 -p 2 



1/2 



cos < 



l-p? 



1/2 



sine 



XlX 2 



(34) 



In fig. 3 we plot the surfaces of constant bi- 
spheroidal coordinate in the x — z plane. In either 
the p ~ 1 or p <C 1 limits, the surfaces of con- 
stant (xi, X2) are nearly in the z and R directions 
over most of the star, as one would expect for a 
local plane wave propagating in the z or R direc- 
tion. (Here R is the cylindrical radius.) However, 



11 One must use care evaluating k r for inertial modes near 
the surface of the star, since it can vary strongly with the 
angle 8. Qualitatively, this strong variation occurs because 
one is imposing a spherical boundary condition on waves 
with inherent cylindrical symmetry. 

12 If the density profile near the surface is a power law with 
depth, one can separate variables in the bi-sphcroidal co- 
ordinates introduced below. These more rigorous solutions 
close to the surface agree with the cutoff behavior described 
here for the density. Although one could, in principle, 
match the interior WKB solution to the surface solution, 
the cutoff for the density gives an adequate approximation 
for the problem at hand. 



the coordinate lines near the point r ~ and 
cos6> = ±\p\ on the surface of the star vary quite 
rapidly with respect to r and 9. The result is that 
the WKB wavenumber becomes quite large near 
these singular points. (We give a detailed math- 
ematical discussion in appendix B.) As one can 
see from fig. 3, the coordinate lines come closer to 
the surface near the singular points, implying the 
upper turning point is much closer to the surface 
near the equator (for small p) than the poles. As 
a result, the wave amplitudes will be much larger 
near the equator, as we will now show. 

Our approximate solution for the interior of the 
star is to ignore terms of order (|fc|£f) -2 , so that 
the differential equation becomes 



0. 



(35) 



Changing to bi-spheroidal coordinates in eq. 35 
gives separable differential equations (see Bryan 
(1889) and Lindblom & Ipser (1999) for details) 
. Define the solution f{x\,x 2 ) — fi(xi)f 2 (x 2 ) 
where both f\ and f 2 satisfy 



d_ 

dx 



(I-* 2 ) 



2^/ 



dx 



m 



f 



((86) 



for separation constant k 2 . This equation has 
the Legendre function solutions k 2 = n(n + 1), 
/i(a;i) = Pnm(xi) and f 2 {x 2 ) = P nm {x 2 ). 
The resulting solution for tjj is then 



ip(x,t) = ip 



1/2 



Pnm {% 1 ) Pnm 2 ) 6 



im<p 



Note the important fact that this solution is valid 
for an arbitrary density profile p, so long as one is 
safely in the short wavelength limit. This is true 
even when p is not a separable function of x\ and 
x 2 , as is generally the case in the interior since 
p(r) = p ^l-(x 2 -p 2 )(p 2 -x 2 )/p 2 (l-p 2 )) . 

The r-modes do not have short wavelengths and 
hence cannot be described by the above WKB ap- 
proximation. However, in the leading order ap- 
proximation of a spherical background star with 
no buoyancy, the r-mode solutions are given by 
(Bryan 1889) 

1p{x,t) = ^oP\m\ + l,m( X ^ P \m\+l,m(x 2 )exp(im(f>- 

oc zi?' m ' exp(im<f> — iuit) 
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Fig. 3. — Surfaces of constant bi-spheroidal coordinate for different values of q. The surface of the star is 
the thick solid line on the unit circle. Only the portions of the spheroids inside this circle are relevant. The 
short dashed lines represent surfaces of constant bi-sphcroidal coordinate — < x 2 < l/|<z| while the long 
dashed lines are for the second bi-spheroidal coordinate, which takes on the range l/\q\ < x\ < 1. The level 
surfaces are at the values 0.1, 0.2, 0.9. The surface of the star is described by two coordinate patches, 
x\ = l/\q\ near the equator, and x 2 = l/q near the pole. 



and have frequencies p = — sign(m)/(|m| + 1). 

We now derive the WKB limit for the solution 
in eq. 37. Writing x\. 2 — cos 6*1,2 and substituting 
/ oc sin' m ' 9 exp(z J e dO k) in cq. 36, we find the 
following standing wave solutions: 

f{6) ~ -^— cos (p0 + a ) (39) 
sin ' 9 

where the wavenumber is given by 

p = (n(n + 1) - \m\[\m\ + 1]) 1/2 ~ n (40) 
for n > |m| (the WKB limit) and 

a = —pir/2 for even parity 

a = — (p + l)7r/2 for odd parity modes(41) 

We have chosen to normalize the Legendre polyno- 
mials to unity over An steradians. Note that the 
nodes are spaced evenly in 9\ t2 . This WKB ap- 
proximation to the Legendre equation fails within 
about one wavelength of 9 = 0, n. The sin #1,2 
factor causes an increase in amplitude toward 
sin 01,2 = 0. The collected result is then 

= to cosjpOi + a) cos(p6> 2 + a) 

n2 sin^sin^V 72 



The factor in the denominator p sin 9\ sin 9 2 oc pR 
is just the mass element, and enforces roughly 
equal kinetic energy in between each pair of nodes. 

An approximate dispersion relation is easily 
derived using the eigenfunctions of eq.42. The 
boundary condition is that the compression term 
£ r / H in eq. 28 must remain finite as H — > 0, im- 
plying £ r — > in the low frequency approximation. 
At either surface patch xi — \p\ or |x 2 | = \p\, this 
condition implies 

(l-x 2 )^+mP nm = (43) 

at x = ±\p\- Eq.43 is equivalent to the one given 
by Bryan (1889), as noted by Lindblom & Ipscr 
(1999). The dependence of the frequency and 
wavenumber on the background stellar model, as 
discussed by Lockitch & Friedman (1999), is small 
in the WKB limit. Substituting the WKB expres- 
sions gives 

Tfl 

sin 9 t&n(p9 + a) = . (44) 

P 

In the limit p ^> \m\, the solutions are found by 
inspection to be p9 + a = —kir, for the mode index 
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k. Including the finite m term to first order gives 
the solution 

= sin(/c-7r/p + m/p 2 — ir5/2p) (45) 
~ kir/n (46) 

where S = for even parity modes and 5 = 1 for 
odd parity modes. The m term can be dropped ex- 
cept for the very low frequency, even parity mode 
with frequency fi ~ m/p 2 . All other modes have 
/i oc p^ 1 . We find the approximate formula in 
eq.45 to agree quite well with the exact solutions 
of eq.43 even for p as small as 5. Lastly, we note 
that Lockitch & Friedman (1999) have checked the 
eigenmodes found using bi-spheroidal coordinates 
with those from a numerical code in spherical co- 
ordinates, finding agreement. 

We choose to normalize the eigenfunctions so 
that at unit amplitude (A = 1) all modes have the 
same energy, which we call 2_E un ; t . We can an- 
alytically compute the mode energy in the WKB 
limit where the eigenfunctions are rapidly oscillat- 
ing (n»ffl), with the result 

E = u 2 jd 3 xp\t\ 2 

- ^- \ x % f ,^ ^ 2 . (47) 

This formula agrees well with numerical integra- 
tions. Our normalization convention is that at 
unit amplitude all modes have the same energy 
2-Eunit- We then find the value of the normaliza- 
tion constant 

2tt ( i_^)3/2 Pofi 2 r-* 

Modes with rapid spatial variation (p > 1) or 
larger frequency \x have smaller normalization in 
order for the energy to be the same. As fi — > 1, the 
wave amplitude goes to zero since inertial modes 
do not exist outside this range. 

Before moving on to discuss the nonlinear force 
coefficients, we discuss the normalization integral 
in eq.47. One can easily find the mode energy to 
leading order in \i by setting fi = in eq.47. In this 
limit, the bi-sphcroidal coordinates become x\ ~ 
(1 — R 2 ) 1 / 2 and X2 — 0. In this limit, the integrand 
is constant in z, and varies as (1 — R 2 )~ x / 2 with R, 
which is large near the surface. The kinetic energy 
then converges as (1 — R 2 ) 1 / 2 from the surface. 



4. coupling coefficients 

The lowest order nonlinear interaction couples 
three inertial waves, implying quadratic nonlinear 
terms as in eq.2. The expressions for the nonlin- 
ear force coefficients can be derived cither from an 
action principle (Newcomb 1962; Kumar & Gol- 
drcich 1989; Kumar & Goodman 1996) or directly 
from the equation of motion (Schenk et al. 2002). 
Note that Schenk et.al. have stressed that the 
form of the coupling coefficient is the same for ro- 
tating systems as for nonrotating systems; only 
the explicit expressions for the eigenfunctions and 
background stellar model need be modified. Since 
we are using daughter modes with wavelengths 
much smaller than a stellar radius, we keep only 
the largest 13 term in the coupling coefficient in 
an expansion of (|fc|_ff) _1 . For modes £ 1; £ 2 , £ 3 , 
the dimensionless coupling coefficient 14 is 



^123 




+ Q&Pi-,ij+&{SP2-,ii)- (49) 

Since Sp cx fi 2 , we find that the interaction en- 
ergy, ft-Eunit; scales as the rotational energy of the 
star. A natural unit of energy is then £? U nit = 
0.5Mr 2 O 2 . In these units, A 2 is the mode energy 
in units of 2£ unit = Mr 2 2 . 

In section 4.1 we discuss conservation rules for 
the nonlinear coupling coefficients. Effectively, 
these rules pick out the largest possible coupling 
coefficients. The scalings for k are discussed in 
section 4.2. We confirm a result found in previ- 
ous studies (Wu & Goldreich 2001) that for waves 
which satisfy the conservation rules, the coupling 
coefficients do not become smaller as the daughter 
mode wavenumber is increased; even though each 
individual eigenfunction is highly oscillatory, the 
product is relatively constant. Numerical results 
are presented in section 4.3 and a detailed analytic 
calculation is given in appendix B. 

Inertial waves in an infinite homogeneous, incompressible 
medium have a nonlinear coupling as given here. Including 
terms arising from compressibility or variation of the back- 
ground stellar quantities then gives terms which are small 
in the limit \k\H 2> I- We ignore these terms here for sim- 
plicity, although the r-mode is formally a large lengthscale 
mode. 

see Schenk et.al. for a derivation of eq.2 and the explicit 
form of the dimensionless coupling coefficient re 
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Fig. 4. — Coupling coefficient as a function of parent mode quantum number n\. The other parent mode 
quantum numbers are held fixed at k\ = 1 and mi = 2. The daughter mode quantum numbers n 2 and 
n3 are allowed to vary, but we have fixed k 2 = —2, m 2 = —3, k% = —1, and raj, = 1. The dots show the 
coupling coefficient determined by numerical integration as described in the text. The line gives the analytic 
approximation k 12 3 = n x [i{ 2 — 7r~ 2 nffc^ 2 . 



4.1. energy and momentum conservation 

Consider a parent mode with quantum num- 
bers (ni,&i,mi) and frequency We are 
free to choose daughter mode quantum numbers 
(n 2 , k 2 ,m,2) and (ri3,fc3,m3) in order to find the 
largest coupling coefficient (see e.g. Wu & Gol- 
dreich (2001)). The integrand is highly oscillatory 
unless the phases of the waves match at each point 
in the star. If we expand the standing wave solu- 
tion in cq.42 in terms of travelling waves, a non- 
oscillatory integrand implies momentum conser- 
vation for the three travelling waves. In addition 
to conservation of the m quantum number, due 
to axisymmetry of the background star, we also 
have momentum conservation along the 9 1 and 9 2 
directions. For small [i, the total number of nodes 
along 6*i and 9 2 simplifies to Ni ~ p ~ n and 
N 2 ~ Mp/tt ~ k. The approximate conservation 
laws which lead to large k can then be written 

mi + m 2 + m 3 = angular momentum 
\n 2 — n 3 \<m momentum along 0\ 
\ — \kz\ \ ^ | fci| momentum along 9 2 . (50) 



For small frequency, the 9\ and 9 2 directions lie 
nearly along the R and z directions, so that the 
second and third momentum conservation rules 
correspond to conservation of momentum along R 
and z. In the limit that the daughter modes have 
much smaller wavelengths than the parent mode, 
which will turn out to be the important limit, we 
find the simple result n 2 ~ n 3 and \k 2 \ ~ |fc 3 |; 
momentum conservation implies that the daugh- 
ter modes have momenta of equal magnitude and 
oppositely directed. 

So far, we have used momentum conservation 
to determine three of the six daughter mode quan- 
tum numbers. In order for energy to be efficiently 
transferred between modes, the interaction must 
be as nearly resonant as possible, meaning that 
the detuning is small: 

5uj = u\ + uj 2 + uj?, — 0. (51) 

There are two simple limits of interest. For short 
wavelength daughter modes with n 2 ~ n 3 and 
k 2 ~ fc 3 , one has [i 2 ~ yU 3 ~ —fii/2; the par- 
ent mode interacts with nearly identical daughter 
modes of half the frequency of the parent. The sec- 
ond solution is where m ~ n 2 , and n 3 -c n\,n 2 . 
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In this case we find ri2 ~ n\ + n^, k2 — k\, 
and k% ~ ki(nz/ni) 2 . The frequencies are then 
p 2 ~ -fii and fi 3 ^ Ati(ri3/ni) < Mi- 

4.2. analytic estimates 

Here we give a back of the envelope estimate for 
the coupling coefficient, leaving the more detailed 
calculation for appendix B. We will only consider 
the important limit of short wavelength, nearly 
identical daughter modes with p, 2 — M3 — —fii/2. 
We shall ignore factors of order unity for the 
present, concentrating only on the scalings. As 
k is dimensionless, we set r* = 1 in this section 
for simplicity. 

Incomprcssibility of the waves implies 



1^2 1 ~ | Mi I an d Zi ^ ^2, the final result is then 



E, 



— [ d 3 x pu\k\ z ^l z . (52) 



For the polytrope of index 1 we find p oc r. t — r = z 
near the surface, where z is the distance from the 
surface. Since the WKB envelope of the waves 
rises steeply toward the surface, and the factor 
of p cancels the p _1 from ^zi we nn d that the 
dominant contribution comes above the turning 
point for the parent mode, where 



ipi 



1 



2 -1/2 ' 



(53) 



Here z\ is the turning point depth of the par- 
ent mode. The daughter mode eigenfunction is 
strongly peaked in the 9 direction due to the 
wavenumber 



k 



2z 



(54) 



where 9 is the polar angle in spherical coordinates. 
The displacement for the daughter mode is then 



6 



z- 1 ' 2 [(cos 2 0-^) 2 + 8Mpr V2 (.55) 



For cos9 ~ \pb2 1, kiz — pi since it is well away 
from the singularity for mode 1 at cos 9 = \(J,i\. 
Plugging into eq.52 we find 



.-1/2 f Sl „ f 1 

K ~ p\z 1 'I dz I 

Jz 2 J-l 



d(cos9) 



(cos 2 9 - pi) 2 + 8p 2 z 



.(56) 



The integrand has a width d(cos0) ~ z 1 / 2 and a 
height (/xp) -1 , giving an area (/ip 1 / 2 ) -1 . Using 



Pi 



dz 



,,2 ? 1/ 2 /- Z 1 / 2 

PlZ^ J Z 2 " 



■ P1IJ>1 



(57) 



The detailed calculation in appendix B confirms 
that the coefficient is about unity. 

We now comment on the scalings for the max- 
imum coupling coefficient in cq.57. The maxi- 
mum coupling coefficient is found to be indepen- 
dent of the daughter mode quantum numbers. The 
reason, elucidated by Wu & Goldreich (2001), is 
that one is integrating over the daughter mode ki- 
netic energy pp^2- This quantity is normalized 
to 2£" un it when integrated over the whole star, 

1 12 

and is roughly £7 un it x Zi when integrated over 
< z < z\. Next, the factor pi ~ n\ implies 
shorter wavelength parent modes interact more 
strongly. This factor would appear for coupling 
of local waves in a box. However, the factor p^ 2 
would not appear for local waves in a box; it arises 
from the large peak in the integrand near the sur- 
face. 

One might wonder whether or not the ap- 
proximate conservation laws for colliding WKB 
waves will hold since one is integrating over a 
small region of the star. The dominant contri- 
bution to the integrand comes from a region of 
size z <~ z\ <~ |/Ui|/ni, and the angular size is 
d(cos6>) - z 1 ' 2 - (l^il/m) 1 / 2 . The daughter 
modes have wavelengths 712 or P2 n 2i depending on 
direction, so there are still sufficient oscillations in 
the important region of the star for large enough 
n 2 . 

4.3. numerical calculation 

We compute the integral in eq.49 numerically 
as follows. Choose a point in the star at which to 
evaluate the integrand. Evaluate ip and Sp on the 
vertices of a Cartesian cube about this point. The 
derivatives in eq.29 and 49 can then be taken along 
Cartesian basis vectors 15 , and then appropriate 
sums over indices taken. The resulting scalar in- 
tegrand is independent of the coordinate 4> since 



15 We evaluate vector quantities along Cartesian basis vectors 
to avoid "curvature terms" (Wu & Goldreich 2001) arising 
from differentiating curvilinear basis vectors. Wu and Gol- 
dreich found these terms are quite large, and cancel out 
in the end, so that significant cancellation error can occur. 
We avoid such cancellation error by using Cartesian basis 
vectors. 
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mi +m2 + TO3 = 0, so that only a two-dimensional 
integral over r and 9 remains. We perform this in- 
tegration with second order accuracy, and increase 
the number of grid points until the integral con- 
verges. 

In fig. 4, we show the numerical integrations for 
the coupling coefficient as a function of m, but 
fixed ki and mi. We also fix (^2,^2) and (£3, to 3 ) 
but allow n 2 and n 3 to vary. For a given m we 
see there is a variation in n due to the degree of 
momentum conservation. However, the upper en- 
velope set by the maximum coupling coefficient 
agrees to within ~ 10% of our analytic formula. 

5. damping and driving rates 

We review the driving rate by gravitational ra- 
diation, and derive simple analytic estimates for 
the damping rates of inertial modes. 

5.1. driving rate 

Gravitational radiation reaction is a driving 
force if the phase velocity in the azimuthal direc- 
tion is positive in the inertial frame and negative in 
the rotating frame; otherwise it damps the mode 
(Friedman & Morsink 1998). The driving rate 
falls off extremely rapidly with wavenumber, so 
that only the very lowest modes have an apprecia- 
ble driving rate compared to damping. Lockitch 
& Friedman (1999) have computed these driving 
rates for the inertial modes of a polytropc of index 
1, and identified several low order driven modes. 
However, since the most unstable mode by far is 
the (n, k, m) = (3, 1, 2) r-mode, we can ignore all 
the others to a good approximation. 

The driving rate of the (n,k,m) = (3,1,2) r- 
modc for a polytrope of index 1 with M = 1.4M Q 
and r* = 12A;m is (Lockitch & Friedman 1999) 
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5.2. bulk viscosity damping 



We now compute the damping rate of inertial 
modes by bulk viscosity damping due to the mod- 
ified URCA processes. We take the coefficient of 
bulk viscosity from Sawyer (1989) and Cutler et al. 
(1990). 

The damping rate is 
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The Lagrangian compression is 
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where the second equality is for low frequency 
modes. The bulk viscosity coefficient is 



C = CfidU 2 T 9 V 



(61) 



where 



Q id = (6 x 10 25 g cm- 1 sec- 3 )(10 15 g cm" 3 )- 2 

(62) 



6 x 10 5 g 1 cm 5 sec 3 . 



Plugging in gives 



-Ebuik = CfidT 9 d x 
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We will evaluate this integral for a polytrope of 
index 1. In this case 



rig gm 

P 2 r*p 



(64) 



is a constant so that 

-E b uik = Cf^ T l(%^) 2 J dWl£r| 2 (65) 

In the WKB limit, this integral is logarithmically 
divergent at r = r, and cos 9 = ±|//|. This diver- 
gence implies that equal contributions to the inte- 
grand come per decade of distance from the sur- 
face. Since the true eigenfunctions flatten off one 
wavelength from the surface, we cut off the inte- 
grals at this distance. Plugging everything into the 
integral and approximating slowly varying quan- 
tities by their surface values gives the amplitude 
damping rate 
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where A = 2|/i|(l — /i 2 ) 1 / 2 ^ is roughly the number 
of nodes along the rotation axis. Evaluating this 
expression for a fiducial neutron star with poly- 
trope index n = 1, mass M — 1.4M Q and radius 
r* = 12km we find the numerical value 



Ibulk = 1.7 x lO-^scc- 1 T Q V 2 



In A 

A* 



Note the extremely important fact that this damp- 
ing rate is very weakly dependent on the wave- 
length of the mode! The usual picture of a cascade 
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to small scales does not make sense for damping by 
bulk viscosity. Instead one must carry the energy 
to small frequency. 

For the (n,k,m) — (3,1,2) r-mode, the previ- 
ously calculated value is (Lindblom et al. 1999) 

lbuik{r - mode) = 2.8 x 1(T 12 sec" 1 T 9 6 t$g8) 

The r-mode has a different scaling with ft and 
normalization since the compression is smaller: 
V • £ oc n 2 instead of V • £ oc f/if. 

5.3. shear viscosity 

The shear viscosity for nuclear matter has been 
calculated by Flowers & Itoh (1979), with an an- 
alytic fit by Cutler et al. (1990) of the form 

v s = VsMp/po) 5/4 Ti 2 (69) 

where v aJld = 2000 cm 2 sec- 1 (p /10 15 g cuT 3 ) 5 / 4 . 
The shear viscosity damping is then 

-E shear = T 9 - 2 J d?xpv s w 2 (j£ (jJ )| 2 - i|V • £\ 2 ^j 

Ts 2 v sJidPo uj 2 J d 3 x(p/ Po ) 9/4 k 2 \Z\ 2 (70) 

where we have kept terms of leading order in 
(|fc|f/") _1 , and subscripted brackets denote a sym- 
metrized derivative. Plugging everything in, and 
approximating the density by a power law with 
depth appropriate for a polytrope of index 1, we 
find the damping rate 

IsHear * 0.6*^ (71) 

For our fiducial star this becomes 

v 2 

Ishear = 3.8 X lO" 9 SCC" 1 T g 2 > j^—j . (72) 

The previously computed r-mode shear damp- 
ing rate is (Lockitch & Friedman 1999) 

lshear(r - mode) = 4 x 10~ 9 sec" 1 T 9 " 2 (73) 

which is about a factor of two different from our 
formula. 

As first noted by Bildsten & Ushomirsky 
(2000), the r-mode is damped much more effi- 
ciently by shear in the crust-core boundary layer 



than by shear over the bulk of the stellar inte- 
rior. Levin & Ushomirsky (2001) later corrected 
this damping rate to account for crust with a fi- 
nite shear modulus. The key parameter is the 
fractional velocity jump over the boundary layer, 
called r\. Levin and Ushomirsky found the rate of 
damping to be 

7„ K (r - mode) = 1.5 x lO^scc" 1 rfuJ^Vji) 

with a realistic estimate for the fractional velocity 
jump of r\ ~ 0.1. Inclusion of the finite shear mod- 
ulus of the crust gives much better agreement of 
the r-mode instability curve with the observations 
of LMXB's. 

We have neglected damping of the daughter 
modes by shear in the boundary layer. 

6. r-mode saturation by discrete modes: 
the small driving limit 

A fundamental plot for the r-mode instability 
is given in fig. 5. The r-mode is unstable for spin 
frequencies above the thick dashed lines, where 
Igr = Ishear, Ivbh or Ibuik- The solid lines show 
where driving of the r-mode equals damping of 
daughter modes, indicating marginal stability of 
the energy transfer. For bulk viscosity, the ratio 
of driving to damping is 

C& = 1*- = 8 x 10% 1Z T 9 - 6 (75) 

Ibulk 

while for shear viscosity 

<U = ^ = 1.3xlO% 6 hz T 2 n- 2 . (76) 

^shear 

In these estimates we have used n = 1/6, appro- 
priate for daughter modes with the largest cou- 
pling coefficients, and n denotes the wavenumbcr 
of the daughter mode. Only in the region from the 
Cs = 1 and C& = 1 lines to the r-mode instability 
curve can we possibly have stable energy transfer 
for the three mode system. 

6.1. young neutron stars 

Nascent, rapidly rotating neutron stars cool 
into the region of instability (Owen et al. 1998) at 
fixed spin frequency. For daughter modes mainly 
damped by bulk viscosity, there is a narrow region 
near the instability curve in which energy trans- 
fer for a single triplet of modes would be stable. 
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Fig. 5. — Stability of energy transfer from the r-mode to daughter modes damped by bulk or shear viscosity. 
The heavy dashed lines show where the r-mode is marginally stable by equating j gr to either ^ shear , Ivbl for 
rj = 0.1, or "fbuik', the r-mode is unstable above these lines. The light solid lines show where "f gr equals the 
bulk (£& = 1) or shear (( s = 1) viscosity damping of daughter modes. The three lines ( s — 1 correspond to 
the number of nodes n = 1, 10, 100 from bottom to top. The stability of energy transfer to discrete daughter 
modes is only stable in the region from the r-mode instability curve (dashed lines) to the ( o s = 1 lines, 
i.e., energy transfer to a small number of discrete modes is not stable in the central portion of the r-mode 
instability region. The nearly vertical dot-dashed line is where bulk viscosity damping equals shear viscosity 
damping for large lengthscale (n = 1) mode. 




However, in this region, the damping is relatively 
independent of the daughter mode wavenumber. 
The quality factor of a daughter mode is roughly 



3.4 x 10 5 ^ hz r io 6 (77) 



where we have used the daughter modes with the 
largest coupling coefficients so that ^ 2 ,3 = 1/6. 
We have also set In A <~ 1. We found the coupling 
coefficients are roughly n ~ 27 for a parent mode 
with m — 3 and = 1/3 so that the saturation 
amplitude for a three mode system is given by 



A\{bulk) = 



1 



E(r — mode) 
0.5Mr2f2 2 ~ n 2 Q 2 d 



~ 10 



-14,-6 rpl2 
''khz 1 10 ' 



(78) 



This formula would imply that nascent neutron 
stars cooling into the instability curve after a su- 
pernova will saturate at a very small fraction of 
the rotational energy of the star. 

However, this formula is not applicable for the 
following reason. Since the damping rate of the 
daughter modes is relatively independent of the 
daughter mode wavenumber, all daughter modes 



have essentially the same parametric threshold 
(roughly eq.78) until n becomes large enough that 
shear viscosity becomes comparable to bulk vis- 
cosity (see fig. 5). We estimate this point to be 
at 7ib= s — lO^i^khz- For daughter modes with 
frequency ji ~ 1/6, there are roughly n^_ s ~ 
10 8 Tf f 1 ^ lz daughter modes parametrically excited 
to large amplitude by the parent r-mode, so that 
the discrete limit is not applicable. Thus the r- 
mode instability in young neutron stars is in the 
continuum limit discussed in section 7. 

6.2. LMXB's 

For the LMXB case, neutron stars with tem- 
perature T ~ 3 x 10 8 K spin up until they hit 
the instability curve (Bildstcn 1998; Levin 1999). 
When the instability curve near T 9 = 0.3 is set 
by boundary layer shear viscosity with r) = 0.1 
(Levin & Ushomirsky 2001), we see that if the star 
stays close to the instability curve, one must go to 
daughter modes with 10 — 100 nodes before the 
energy transfer can become stable. As a result, 
10 2 ~ 4 daughter modes will be parametrically ex- 
cited to large amplitude, and the continuum limit 
is more appropriate. 
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If, however, boundary layer shear viscosity is 
not operating for some reason, then the discrete 
mode approximation will be valid near the insta- 
bility curve. The quality factor of the daughter 
modes in this case is 



Q d ~ 5.5 x lOVhzTf n" 



(79) 



giving a saturation amplitude for LMXB's near 
the instability curve to be 



A 2 (shear) = 



E(r — mode) 1 

0.5Mr 2 il 2 ~ k 2 Q 2 
10- 21 (0.33/^ khz ) 2 T- 4 n 4 ,(80) 



which is quite small. 

The conclusion we draw in this section is that, 
for the likely scenario in which either bulk viscosity 
or boundary layer shear viscosity sets the r-modc 
instability curve, many modes will be parametri- 
cally excited to large amplitude, and the contin- 
uum limit discussed in the next section is a better 
approximation. 

7. r-mode saturation in the continuum 
limit 

The saturated r-mode energy was found in sec- 
tion 2.2 to be 
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This solution is valid when there is a clear sepa- 
ration between the inner and outer scales of the 
turbulence (see Fig. 2). The outer scale is given 
by the r-mode itself, while the inner scale is where 
7ni, the characteristic rate for amplitude change 
by nonlinear interactions, is equal to the dissipa- 
tion rate given by 7 s hoar (bulk viscosity is irrel- 
evant for the inner scale; see below). In other 
words, the inner scale is where the Reynolds num- 
ber for that scale becomes order unity. We can es- 
timate the nonlinear interaction rate using cqs.12 
and 13, with the scalings n = n a ~ np ~ n 7 , 
M = Ma ~ A*/3 ~ M7) etc - We find 

7ni ~ -rz - -j- T — OiV 2 ~ 7gr — 82 

where n r and pL r set the scale for the driving region 
at which 7 gr = 7 nl . 



The expression for shear viscosity from eq.72 
can be written 7 s hoar = Jsn 2 in the // 2 <gc 1 limit. 
Equating j n \ to 7 s hcar, we find the inner scale is 
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In the region of the (f2, T) plane where the right 
hand side of Eq. (83) is large, many modes are 
excited with a clear separation between inner and 
outer scales of turbulence (see Fig. 2) . This region 
includes the entire instability window when the 
r-mode is damped by a viscous boundary layer. 
When there is no viscous boundary layer, there is 
a small region close to the instability curve where 
(83) is small and where the discrete limit applies 
instead of the continuum limit. 

For bulk viscosity things work differently. The 
cascade solution will be valid in the region of phase 
space where 7 nl > 7 bu ik, or 
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where we have written 7b u ik = lb^~ 2 - In the re- 
gion of the (f2, T) plane where the right hand side 
of Eq. (84) is large compared to unity, bulk viscos- 
ity is dominant at the outer scale and no cascade 
solution exists. When the right hand side of Eq. 
(84) is small compared to unity, then a cascade 
can form, but the bulk viscosity is irrelevant for 
setting the inner scale of the cascade. We note 
that the boundary (84) of the bulk viscosity dom- 
inated regime approximately coincides with the 
curve Cfi = 1 in Fig. 5. A newly born neutron star 
will very rapidly move from the instability curve 
to the C& = 1 curve at which point the cascade can 
form. 

Finally we note 16 that the weak turbulence ap- 
proximation which underlies the derivation of Eq. 
(12) eventually breaks down as one goes to small 
scales. The weak turbulence approximation re- 
quires that the nonlinear energy transfer timescale 
l/7ni ~ 20 s J / khz rt_1 M 3 ^ 2 be much longer than the 



mode period 



5 x 10 4 s^i 1 ^ kh 1 z , which breaks 



3 We thank P. Goldreich for bringing this to our attention. 
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down in the regime 

n^ -5/2 >4x lO 4 ^. (85) 

Thus our turbulent cascade solution of the equa- 
tions of motion (12) will likely be replaced by some 
form of strong turbulence at sufficiently small 
scales. However, this should not affect our pre- 
diction of the r-mode's saturation amplitude, as 
the regime (85) in phase space where the approx- 
imation breaks down is well separated from the 
driving regime fj, ~ n ~ 1. 

8. comparison with previous work 

There have been three distinct alternative 
nonlinear mechanisms proposed to saturate the 
growth of the r-mode: (1) For large amplitude 
pulsations, the Fermi energies of the electron- 
proton-neutron gas become significantly shifted, 
and the kinetic energy can be rapidly converted to 
both heat and neutrinos by nonlinear bulk viscos- 
ity (Reisenegger 2001). (2) The amplitude grows 
so large (E ~ E Iotation ~ GM 2 /r,) that strong 
shocks occur, rapidly thermalizing the kinetic en- 
ergy (Lindblom et al. 2001, 2002); (3) In neutron 
stars with a crust, a turbulent boundary layer 
forms at the crust-core interface (Wu et al. 2001). 
We now discuss each in a bit more detail. 

Since the Fermi energy of the electron and 
neutron have a different dependence on density, 
the Fermi surfaces are shifted out of beta equi- 
librium when matter is compressed. The scal- 
ing of the resulting neutrino emission rate de- 
pends on the ratio of chemical potential imbal- 
ance, Ap,„ = ii n — Hp — jJL e , to temperature, which 
is roughly (Reisenegger 1995) 

T 3 T p V 1 

where Ef.c is the Fermi energy of the electron. 
When this ratio is large, 5/8 of the resulting dis- 
sipation heats the star and 3/8 goes into neu- 
trinos. The rate of such dissipation scales as 

( A| f 33 T ) (X (^p/p) S times the neutrino emissiv- 
ity of uncompressed matter. The mode damping 
rate is extremely sensitive to the compression, and 
can saturate the growth of the r-mode for suffi- 
ciently large amplitude. Reisenegger (2001) has 
done a detailed calculation, finding that the satu- 
ration energy is comparable to the stellar rotation 



energy. This interesting idea gives a larger (less 
constraining) saturation amplitude compared to 
the value found in this paper. 

Next, Lindblom et al. (2001, 2002) have per- 
formed state-of-the-art 3D Newtonian hydrody- 
namics simulations including a prescription for the 
radiation reaction force. The only damping mech- 
anism included in the code is numerical viscos- 
ity, and of order 128 3 points were used. They 
were able to follow the linear growth of the r- 
mode, all the way into the nonlinear regime where 
E > E rot . Shocks then formed near the surface of 
the star, rapidly thermalizing the kinetic energy 
of the mode. 

Since the growth of the instability is so slow 
compared to the dynamical time in the star, they 
found it necessary to artificially increase the ra- 
diation reaction force by a factor of <~ 4500. A 
natural question is how the mode would saturate 
if the correct, physical value of the driving force 
was used. The following physical example is useful 
to consider. Imagine water waves being driven by 
wind moving at 1 cm sec -1 , a whisper of a breeze, 
as compared to 4500 cm sec -1 , a hurricane. For 
small amplitude water waves, four wave interac- 
tions can transport energy to small scales, satu- 
rating the growth of the waves. In a hurricane, 
the wave growth time is so short that waves grow 
to large amplitudes and break. Since Lindblom et 
al. have not addressed how saturation might oc- 
cur for physical values of either driving or damping 
of the waves, the relevance of their simulations to 
saturation of the r-mode instability is not clear. 

One comparison which can be made is to 
use our formula in eq.20 to estimate the sat- 
uration amplitude seen in Lindblom et al.'s 
simulations when -f gr — > 45007 gr . We find 
£ simu i ation /(0.5Mr2ft 2 ) = Q1 ( ae / .i) (4500 7gr /ft) ~ 
5 x 10 -3 (a e /0.1)^ hz . This result can be trans- 
lated into Lindblom et al.'s notation by using 
A\ ~ 0.5Ja s 2 imulation = 0.008a| mu , at ; on with the 

result a s i mu i a tion ^ 0.7 (a e /0.1) 1/2 i^hf - This com- 
parison shows that if one attempted to extrapolate 
the saturation amplitude over more than three 
decades in driving force, that the saturation am- 
plitude by mode coupling would indeed be of order 
unity, in their notation. This does not, of course, 
explain the saturation amplitude seen in Lindblom 
et al.'s simulations, which they explain is due to 
strong shocks near the stellar surface. However, 
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the results of this paper show that their claims of 
(1) saturation energy of order the rotation energy, 
and (2) strong shocks as the saturation mecha- 
nism, arc not supported. They are an artifact 
of the unphysically large value for the radiation 
reaction force. 

We comment further on the ability of a numer- 
ical simulation to accurately reproduce the cas- 
cade of energy to small scales as derived in this 
paper. In simulations with 128 3 points, only a cer- 
tain number of modes exist because of the finite 
resolution. Since the detuning is a rapidly decreas- 
ing function of wavenumber, secular energy trans- 
fer by nearly resonant interactions becomes more 
important as the number of grid points increases. 
For instance, daughter modes with half the fre- 
quency of the r-mode have quantum numbers in 
the ratio k/n ~ 1/20, so that one needs nearly 
20 times as many nodes in the cylindrical radius 
as in the z direction in order to find parametri- 
cally excited daughter modes. We estimate that 
only a few of these might have been accurately 
modeled by Lindblom et al's simulations. In time 
evolutions of the mode amplitude equations with 
a small number of low order, nonresonant modes 
(Morsink 2001; Arras 2001), large saturation am- 
plitudes were found as compared to the results in 
this paper. The reason, as can be clearly seen 
in eq.5 for the parametric threshold, is that the 
r-mode cannot easily excite daughter pairs with 
large detuning. However, going to higher order 
modes with much smaller detuning can give a satu- 
ration amplitude orders of magnitude smaller than 
for arbitrary, low order modes. 

Lindblom et al. specifically commented that 
three-mode coupling is not the saturation mech- 
anism in section H of their paper. Their claim 
was based on the lack of power observed in cer- 
tain modes besides the r-mode during their sim- 
ulation. However, they focused on interactions 
which couple the r-mode twice to a third mode. 
As they themselves comment at the end of section 
H, they have not included parametric excitation of 
daughter modes in their constraints. As discussed 
in section 2.1, couplings of the type discussed by 
Lindblom et al. (2002) are far less important than 
parametric couplings, because (1) they are down 
by a factor of parent mode amplitude, which is 
small, and (2) only a relatively small region of 
phase space couples well with the r-mode by non- 



parametric couplings. Hence, Lindblom et al.'s 
constraints are not useful since they constrain an 
unimportant process. 

The inability of simulations to include very high 
order modes presumably also explains the results 
of the fully relativistic simulations of Font & Ster- 
gioulas (2001), in which an r-mode with order 
unity amplitude was observed not to lose any en- 
ergy to other modes over several dynamical times. 
More recent simulations by Gressman et al. (2002) 
show that for slightly larger initial amplitudes, the 
r-mode decays rapidly into a differentially rotat- 
ing configuration without shocks forming. These 
results are not inconsistent with our analyses, but 
our results indicate that the r-mode never reaches 
the regime of rapid nonlinear decay seen by Gress- 
man et al. (2002). 

Next we discuss the turbulent boundary layer 
mechanism of Wu et al. (2001) which operates in 
neutron stars with a crust. Energy dissipation 
by turbulent drag scales as A\, leading to sat- 
uration of the mode. The attractiveness of this 
idea is that the turbulent drag force is well under- 
stood in magnitude and scaling both from numer- 
ical estimates as well as laboratory experiments. 
These authors considered the effect of such en- 
ergy dissipation on the crust and thermal history 
of the star, and go on to discuss the observable 
spin frequency of the star after it exits the r-mode 
instability region. For a realistic fractional ve- 
locity jump across the crust-core boundary layer, 
r\ ~ 0.1, they found the r-mode saturated at a 
value E/(0.5Mr*n 2 ) ~ 0.2i/£g z , which is larger 
(less constraining) than the value found here, both 
in normalization, and in the dependence on v\^ z . 
Furthermore, their mechanism does not operate in 
completely fluid stars without a crust, which is the 
case for hot young neutron stars. 

Lastly, we mention that this paper is a compan- 
ion paper to that of Morsink (2002), which dis- 
cusses the nonlinear coupling among r-modes in a 
star for which buoyancy forces are dominant over 
Coriolis forces. Morsink found that, because the r- 
mode frequency decreases with m, interactions do 
not become more resonant as the daughter mode 
m increase. As a result, energy transfer among 
three r-modes is not likely to produce a saturation 
value as low as in this paper. 

We conclude that nonlinear mode coupling to 
inertial modes provides the most stringent con- 
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straints on the r-mode amplitude at this time. 

9. spin evolution of neutron stars 

The spindown torque exerted on the neutron 
star by gravitational radiation is roughly 
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Since the spindown rate decreases strongly with 
spin frequency most of the time is spent at lower 
frequencies. 

The spindown time becomes >10 4 yr at the low- 
est rotation rates inside the instability curve, while 
it is of order a few days for stars rotating near 
breakup. This is of interest for certain gamma-ray 
burst models, such as the "supranova" model (Vi- 
etri & Stella 1999) in which core collapse leads to 
ejection of the stellar envelope, as well as a rapidly 
rotating neutron star which is above the maximum 
mass for a nonrotating star. Angular momentum 
transport can then slow the neutron star down, 
leading to collapse to a black hole and generation 
of a powerful gamma-ray burst. Our results imply 
that the gamma-ray burst should occur within of 
order a week after the supernova explosion. 

Next we turn our attention to neutron stars 
in LMXB's. The ratio of spindown torque, due 
to radiation reaction, to accretion torque r acc ~ 
M(GMr r ) 1/2 is roughly 



Tacc 



2 x 10 42 (a e /0.lKg z 
10 34 (M/10- 8 M© yr- 1 ) 
10- 8 Af Q yr" 1 



M 



For accretion rates smaller than the Eddington 
rate M ~ lO -8 M Q yr -1 , and spin frequen- 
cies above j/ spin = 200 Hz, the radiation reaction 
torque is larger than the accretion torque and can 
halt the further spinup of the neutron star. If 



the neutron star viscosity is dominated by nor- 
mal matter, then the star enters into a limit cy- 
cle of spinup by accretion and spindown by the 
r-mode, as discussed by Levin (1999). [The alter- 
native equilibrium scenario is discussed in Sec. 10 
below.] Since the r-mode is only likely to be unsta- 



ble for 



> 300 Hz, the r-mode can halt spinup 



inside the region of instability. 

The observable spin frequency is determined by 
where the star exits the region of r-mode instabil- 
ity, if no other process spins the star down fur- 
ther. The exact spin frequency at which the star 
exits the region of r-mode instability depends on 
the evolution of both the spin frequency and the 
stellar temperature (Levin 1999; Owen et al. 1998; 
Wu et al. 2001). We can estimate this terminal fre- 
quency (Wu et al. 2001) by equating the neutrino 
cooling luminosity, L u = 7.4 x 10 39 T|ergsec _1 , 
with the rate of stellar heating due to the r-mode. 
If we approximate that all the energy input to the 
r-mode by radiation reaction is damped away as 
heat, the rate of heating of the star is just given by 
-Eheat = 2j gr E, where E is the saturation energy 
found in eq.81. Equating heating and cooling, we 
find the equilibrium temperature as a function of 
spin frequency, given by 



10 9 K 
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The crystallization temperature of the crust is 
(Wu et al. 2001) T mo it - (5 - 10) x 10 9 K so the 
heating by the r-mode cannot prevent the crust 
from forming when ^ sp i n <C 10 3 Hz. If the instabil- 
ity curve is set by boundary layer shear viscosity 
(igr — Ivbi), the intersection of the equilibrium 
spin down curve with the r-mode instability curve 
is given by the terminal frequency 

/ n \ 0.02 / „ n. 0.28 

^ terminal ~ 250 Hz (^J (^-J (91) 
with a core temperature of roughly 

/ n \ 0.15 / v . \ 0.46 

r — = 6 x 10 K (or) (3355) < 92 » 

Note that the observable spin frequency is very 
insensitive to the saturation parameter 
well as to the fractional velocity jump r\. The 
spin frequency found in eq.91 is comparable to 
the lower end of the observed LMXB's, consis- 
tent with a limit cycle (Levin 1999) of spin-up 
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by accretion and spin-down by the r-mode. The 
timescale to exit the instability curve is roughly 
t S pindow„ 2000 yr (330/250) 11 (ae/0.1)- 1 ~ 
4 x 10 4 yr (a e /0.1) _1 . This spindown timescale 
is very sensitive to the position of the instability 
curve. 

For young neutron stars with strong magnetic 
fields, the spindown torque from magnetic dipole 
radiation is comparable to that from gravitational 
radiation. Equating the magnetic dipole spin- 
down timescale (Shapiro & Tcukolsky 1983) i m d — 
30yri?f 2 2 ^ h 2 z to i sp i n down, we find that gravita- 
tional radiation reaction dominates for frequencies 
above v spin ~ 400 Hz (0.1/ae) 1 / 9 B^ 9 , where B 12 
is the surface dipole field in units of 10 12 G. Hence 
for typical pulsars with magnetic fields ~ 10 12 G, 
the spindown torque is dominated by the r-mode 
only for fairly large spin frequencies. 

10. Detectability of gravitational waves 

We now discuss the prospect of detecting grav- 
itational waves from r-modes, based on the satu- 
ration amplitude (21). We consider three different 
scenarios: (i) newly born neutron stars where an 
optically observed extra-Galactic supernova pro- 
vides the sky location for the gravitational wave 
search; (ii) LMXB's in the spinup-spindown limit 
cycle first discussed by Levin (1999); and (iii) 
LMXB's in spin and thermal equilibrium. 

For newly born neutron stars, Brady & Creighton 
(2000) (BC) discuss the detection likelihood by 
LIGO assuming a large saturation amplitude. 
They parameterize the saturation amplitude in 
terms of a parameter n in their Eq. (7.3), which 
we denote by Kb c - Our result (21) gives Kb c ~ 
1.2 x 10 _3 fjj hz , while BC took n bc = 1. In the 
first year of spindown, ^khz decreases from <~ 1 
to <~ 0.66 [cf. Eq. (88) above], and thus the grav- 
itational wave strain amplitude will be a factor 
y/Kbc ~ 2.8 x 10~ 2 smaller than that considered by 
BC. The distance to which the source can be seen 
by enhanced LIGO detectors, for fixed integration 
time (see below) , is correspondingly reduced from 
BC's estimate of ~ 8Mpc to <~ 200 kpc, almost 
inside the Galaxy. Since the galactic supernova 
rate is roughly once per 50 — 100 yrs, the proba- 
bility that LIGO will detect young neutron stars 
radiating due to r-modes is small. 

We now discuss why we can treat the integra- 



tion time as fixed. The matched filtering signal to 
noise ratio S/N for gravitational waves, when av- 
eraged over source orientations and polarizations, 
depends only on the energy per unit frequency 
dE/df of the waves (Flanagan & Hughes 1998): 
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Here D is the distance to the source and Sh(f) is 
the detector noise spectrum. For waves of fixed 
azimuthal quantum number m, using the replace- 
ment dE = 2-irfdJ/m yields (Blandford 1984; 
Lindblom & Owen 2002b) 
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(94) 



where J is the ^-component of angular momen- 
tum. As noted by Lindblom & Owen (2002b), the 
expression (94) is independent of how quickly the 
star looses angular momentum, and hence of the 
saturation amplitude. Thus, a priori one would 
not expect our low saturation amplitude (21) to 
affect very much the detectability of the signal. 
The problem however is that it is not possible 
to integrate long enough to accumulate the total 
signal-to- noise ratio (94). 

Using the stellar model discussed before Eq. 
(58), the relation / = 4^ sp i n /3 between gravita- 
tional wave frequency / and spin frequency, the 
broadband LIGO-II noise curve 17 from Gustafson 
et al. (1999), and neglecting the spin dependence 
of the moment of incrtial of the star, we can eval- 
uate (94) for a spindown from an initial spin fre- 
quency ^khz,i in kHz to a final spin frequency fcfchz.f • 
The result is 
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(95) 



The complete spindown from say 1 kHz to ~ 
250Hz [cf. Eq. (91) above] gives S/N ~ 21 
at lOMpc. The first year of spindown from 
1000 Hz to 650 Hz [cf. Eq. (88) above] gives in- 
stead S/N ~ 6.2, which is not too much smaller. 

However, the need to perform a search over 
spindown parameters in practice limits the inte- 
gration time to ~ 10 6 seconds. BC analyzed the 



17 In the relevant frequency range / > 500 Hz, this noise 
curve is approximately given by fSh(f) ~ 1.7 X 
10- 44 (//1000 Hz) 3 . 
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performance of a "stack slide" search method, in- 
volving demodulating the signal for many different 
choices of spindown parameters, dividing the de- 
modulated data into several chunks or "stacks" , 
computing the power spectrum of each stack and 
adding the power spectra. The threshold value p t h 
of S/N for this method, assuming 1% false alarm 
probability, is approximately given by solving the 
equation 

T(N S , N s + p 2 th /4)/T(N s ) = 0.01/ (N b N p ). (96) 

where T(a, b) is the incomplete gamma function, 
N s is the number of stacks, AT b is the number of 
frequency bins per stack, and N p is the number of 
points in the space of spindown parameters. BC's 
estimate that r-modes are detectable out to 8 Mpc 
was based on assuming a Teraflop of computing 
power, which implied an optimum detection strat- 
egy of ~ 8 stacks of ~ 10 5 s duration each, inte- 
grating from 

^spin — 200 Hz to Z^spin ^ 186 Hz. 

We can modify the BC analysis for our turbu- 
lent cascade scenario as follows. Optimum sen- 
sitivity is achieved late in the spindown, so we 
assume that Z4 p i n ~ 650 Hz, corresponding to 1 
year after the start of the spindown if the initial 
spin frequency is 1 kHz. We take the parameter 
values ^i max = 0.3, / max = (4/3)650Hz (instead 
of 200 Hz as in BC), and a spindown timcscalc 
Tmin = 1 yr, which from Eq. (88) is appropriate 
after 1 year of spindown. Maximizing over the 
number of stacks and stack durations as in BC 
gives that the optimum detection strategy for a 
Teraflop of computing power is to use ~ 10 stacks 
of duration ~ 3 x 10 s each. The corresponding 
number of parameter space points is N p ~ 6 x 10 8 , 
from Eq. (2.20) of BC, which gives from Eq. (96) 
a threshold value of pth ~ 15. Combining this 
with Eq. (95), and noting that at ^ sp i n <~ 650Hz 
an integration time of 3 x 10 5 s corresponds to 
Ai^jpin <~ 0.5(a e /0.1) Hz gives that the source 
would be detectable to <~ 200(a e /0.1) 1 / 2 kpc, con- 
sistent with our earlier estimate 18 . 

This conclusion, however, is based on the as- 
sumption of using the stack-slide search method. 
It is conceivable that an alternative signal process- 
ing strategy (and increased computational power) 

The main reason for the loss of sensitivity compared to BC 
is the reduction in A^ sp i n from ~ 14 Hz to 0.5 Hz; the star 
is spinning down more slowly. 



might enable one to integrate for longer periods 
and achieve a sensitivity closer to the original BC 
estimate. 

We mention in passing another possible dif- 
ficulty in searching for the signal from r-modes 
when a turbulent cascade is present. This diffi- 
culty is that the phase of the r-mode will wander 
randomly in time due to the interaction with the 
turbulent cascade, on some timescale t c . The peak 
in the Fourier transform of the demodulated data 
stream will correspondingly be smeared out over 
a frequency interval of width ~ l/t c , which will 
be over several frequency bins if the stack size is 
larger than t c . 

The phase coherence timescale for a typical 
mode in the cascade will be of order 



or smaller, where 7 n i is the the nonlinear energy 
transfer rate (82) (Zakharov et al. 1992). For the 
r-mode this is only ~ 200 s at v sp i n ~ 700Hz. How- 
ever, one might expect the coherence time of the 
r-mode to be somewhat longer than the estimate 
(97), since the r-mode is being pumped coher- 
ently and is loosing energy by interacting simul- 
taneously with a large number of different modes. 
Unless the phase coherence time for the r-mode is 
10 2 — 10 3 times larger than the estimate (97), the 
sensitivity of the search will be reduced. Again, it 
may be possible to modify the data analysis pro- 
cedure to compensate for the phase wandering 19 . 

Next, we consider the detectability of r-modes 
in LMXB's in the spin up/spin down limit cy- 
cle. The signal from the spin down phase is es- 
sentially the same as for newborn neutron stars, 
except that they will typically be seen at a low 
frequency where most of the spindown time is 
spent. At lower spin frequencies the search over 
spin-down parameters becomes significantly eas- 
ier, since the spindown timescale is longer. The 
formula (4.3) in BC for computational power, with 
/max = fspin and T min given by Eq. (88), shows 
that for ^ S pin ^5 400 Hz integration times as long 
as 10 7 s can be achieved with 1 Teraflop of com- 
puting power. Combining Eqs. (88) and (95) gives 

'The method suggested by BC to compensate for phase wan- 
dering requires a stack size shorter than t c and a compu- 
tational power that scales as 3 T / tc , where T is the total 
integration time. In practice this limits T to < (20 — 30)t c . 
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that the signal to noise ratio for a 10 7 s integration 
starting at ^ sp i n is 




In the regime j/ sp i n < 400Hz, the signal-to-noise 
threshold from the BC method is p t h ~ 10, within 
a factor of <~ 2, giving that the signal should be 
visible to a distance 

/ ;/ . \ 9/2 / n x 1/2 

for !4pin < 400 Hz. 

However, as noted by Levin (1999), the chance 
of observing a particular source emitting gravita- 
tional waves is proportional to the relative length 
of time spent in the spindown phase of the limit cy- 
cle. Using our saturation amplitude we find a duty 
cycle ~ 10~ 3 (0.1/a e ), implying that one would 
need of order 10 3 (a e /0.1) LMXB's within the dis- 
tance (99) in order to overcome the small duty cy- 
cle. Nevertheless, as argued by Heyl (2002), there 
may be enough Galactic LMXB's that some will be 
seen in the spin-down phase by enhanced LIGO, 
especially if a e is smaller than 0.1. 

We note that for LMXB's, the phase wander- 
ing of the r-mode due to the turbulent cascade 
is less of a problem, since the nonlinear energy 
transfer timcscale (82) increases rapidly as f sp i n 
decreases. There is in addition a phase wandering 
due to fluctuations in the accretion torque, but 
this occurs over much longer timescales and can 
be dealt with in the manner suggested by BC. 

The third possibility we consider is when the 
viscosity of an accreting star is independent of 
temperature or is an increasing function of tem- 
perature. In such a case the star can achieve 
an equilibrium state where the accretion spin-up 
torque is stably balanced by the radiation reaction 
torque due to the r-mode, and r-mode heating is 
balanced by neutrino cooling (Levin 1999). Such 
equilibria have been found for stars with hyperon 
cores (Wagoner 2002) and for strange stars (An- 
dersson, Jones & Kokkotas 2001). [However, the 
central densities of neutron stars are sufficiently 
uncertain that hyperon cores may or may not ex- 
ist.] In these scenarios, the equilibrium r-mode 
amplitude is not set by the turbulent cascade con- 
sidered here, but instead by the equilibrium con- 
ditions. The gravitational wave signal is weaker 



than the limit cycle case considered above, and 
its strength can be inferred from the X-ray flux 
(Bildsten 1998). 

11. conclusions 

In this paper, we have accomplished several ob- 
jectives, which can be divided into stellar oscilla- 
tion theory, and phenomenology of neutron star 
spin evolution. 

We have, for the first time, presented a WKB 
theory of global stellar inertial modes (section 3 
and appendix B), including both the rapidly vary- 
ing phase and amplitude which rises quickly to- 
ward the surface. Both the eigenmodes and pulsa- 
tion frequencies take on a very simple form, which 
was never clearly elucidated in previous calcula- 
tions [e.g., Lindblom & Ipser (1999); Lockitch & 
Friedman (1999).] Appendix B gives a detailed 
mathematical treatment of inertial waves. We 
have estimated when the affects of buoyancy be- 
come important. The damping rates by bulk and 
shear viscosity appropriate for neutron stars have 
been derived in section 5, and reduced to simple, 
accurate formulae giving the scalings with stellar 
parameters and mode quantum numbers. 

Next, we have given a complete review of satu- 
ration of an overstable mode in the two different 
limits of strong and weak driving force. The lit- 
erature for the weak driving limit, familiar from 
studies of main sequence or white dwarf pulsators 
(Dziembowski & Krolikowska 1985; Wu & Goldre- 
ich 2001), is reviewed in detail, and shown not 
to apply to most physical situations in which the 
r-mode instability operates. The strong driving 
limit, in which a turbulent cascade forms, has 
never been applied to stellar oscillations to our 
knowledge. Therefore, the weak turbulence meth- 
ods in this paper may find application for ampli- 
tude saturation in stars with a driving force strong 
enough to parametrically excite many modes, but 
weak enough that shocks do not form. In the 
strong driving limit, we find that the r-mode sat- 
urates at an energy E/E rotation < 7 gr /f2. 

The consequences of these calculations for neu- 
tron star spin evolution are as follows. First, 
the time scale for a rapidly rotating (10 3 Hz) 
young neutron star to spin down to a frequency 
Jospin by gravitational radiation spindown torque is 
ispindown ~ 2xl0 3 yr K/O.l^Vspin/SSOHz)- 11 . 
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Hence, the r-mode can be responsible for an ini- 
tial, rapid spindown, but magnetic dipole ra- 
diation will dominate for spin frequencies be- 
low z/ spin ~ 400Hz(0.1/a e ) 1/9 -B^ 9 . Second, 
in spite of the small saturation amplitude, wc 
find that the gravitational radiation spin-down 
torque is still sufficiently large to halt the spin- 
up by accretion for neutron stars in LMXB's. 
Hence, our calculation confirms the validity of 
this assumption by previous investigators. If 
the viscosity is dominated by normal matter, 
the star will enter into the limit cycle in spin 
and temperature discussed by Levin (1999). Fi- 
nally, we estimate that newly born neutron stars 
will be visible to - 200(a e /0.1) 1 / 2 kpc with en- 
hanced LIGO interferometers, and LMXB's in 
the spin down phase of the Levin limit cycle 
out to - 90 (^ S pi n /330 Hz) 9 / 2 (a e /0.1) 1 / 2 kpc for 
i/ spi n<400Hz. 

It is a pleasure to acknowledge many useful 
conversations on stellar oscillations with Yanqin 
Wu. P.A. would also like to thank Chris Matzner, 
Chris Thompson, and Maxim Lyutikov. Part of 
this work was completed when two of us (P.A. 
and S.M) were at the Institute for Theoretical 
Physics at the Spin and Magnetism in Young Neu- 
tron Stars workshop. Wc thank Lars Bildsten for 
his gracious hospitality and for a number of use- 
ful conversations, and John Friedman, Curt Cutler 
and Peter Goldrcich for useful comments on the 
manuscript. P.A. is supported by an NSERC fel- 
lowship. This work was supported in part by NSF 
grants PHY-9900672 and PHY-0084729 at Cornell 
University. E.E.F. was supported by NSF grants 
PHY-9722189 and PHY-0140209 and by the Al- 
fred P. Sloan Foundation. S.M. received support 
from NSERC. 



20 



A. cascade solution 



In this appendix we find scale-free solutions to the kinetic equation 13. The material in this appendix 
can be found in texts on weak turbulence theory, e.g. Zakharov et al. (1992), section 3.3. For convenience 
we use dimensionless units with _E un it = 1 and fi = 1. 

We approximate the mode frequency as uj = 2/i ~ 2i:k/n = s2i:\k\/n, where s = ±1 is the sign of 
the frequency For notational simplicity, we will use both k and /x as positive numbers for the following 
derivation, using s to take into account the sign. The resulting expression will then be written in a form 
valid for either positive or negative H- The sum over modes is then given by 
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Here we used fi — nk/n instead of k since k has an implicit scaling with n. The coupling coefficients are 
assumed to be negligible if momentum conservation is not satisfied, so that 
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\K a i3j\ 2 d(s a in a + spin i3 + s 7 in 7 )(5(s Q2 fca + sp 2 kp + s 72 fc 7 )(5(m ct + mp + m 7 ). (A2) 



The signs s q2 etc. allow for waves moving in either direction (see appendix B). We use the fact that 
the coupling coefficient is approximately independent of both m and the sign of fi, and is separately scale 
invariant in n and fi to say 
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In appendix B we show that u — 1 and v = — 2. Since the mode frequencies and coupling coefficients are 
approximately independent of the m quantum numbers, the m dependence can be integrated over giving the 
function 



dm a I dmp I dm^ S(m a + mp + m 7 ). 

-n a J — np -J —TLry 

This function is symmetric, and is A ~ 4n a np in either the limit n a <C np ~ n 7 or np <C n a ~ n 7 . 
Plugging these definitions into eq.13, and summing over the frequency signs we find 
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We attempt to find scale-free inertial range solutions of the form 

N = N n- p ii- q . 



(A5) 



(A6) 



After inserting this power law dependence on n and \i into the eq.A5, one can make a change of coordinates 
in the second and third sets of terms in the last parenthesis, in order to make them have the same form as 
the first term, up to a scaling factor. For instance, in the second set of terms let np = n^/np', \ip = /j^/hP', 
n 7 = n^'Ha/up', and ^i 7 = n-y'Ha/ Hp' ■ After simplifying and collecting terms we find 

J dm a I a = ttN 2 (2/tt) 2 J dnp npdn^ n^d^p dfi 7 A\R a p 7 \ 2 /j, a upHjS(Sn)d(dk)S(fi a — Hp ~ Hi) 
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where po = u + 3 and g = v + 5/2. We can make the scaling of eq.A7 with n a and [i a explicit by defining 

x f3 = n f3 /n a , x 1 = n 7 /n Q , y (i = /ip/ '/j, a , y 1 = ^/fi a , n a(il = n^ v a f af}l , and A = A/n 2 a . We find 



dm a I a = 4N*n 2 J^-ri- 2 vl {q °- q) - 2 I(p,q) (A8) 

where l(p, q) is the dimensionless integral 

2(p,q) = J dx^x^dx^c^d^A/ 2 ^^ 

x (l - xf p - po) yl {q - qo)+1 - x z(p-Po)y^-io)+^ . (A9) 

Stationary solutions to eq. A7 are now easily found by using the delta functions to force either the first 
or second parenthesis to zero (Zakharov et al. 1992). For instance, using the frequency delta function to set 
the first parenthesis to zero gives N oc fj, , so that all modes have the same energy, i.e., thermodynamic 
equilibrium. Using the momentum delta functions to set the first parenthesis to zero gives thermodynamic 
equilibria with respect to a moving reference frame. Setting the second parenthesis to zero using the momen- 
tum delta functions gives solutions supporting a constant momentum flux. We are interested in solutions 
which support an energy flux. The frequency delta function can be used to set the second parenthesis to 
zero if we let p = po = u + 3 = 4 and q = qo = v + b/2=\/2. We now proceed to show that this solution 
corresponds to a flux of energy to small frequency and large wavenumber. 

In the inertial range, driving and damping are negligible, and the conserved energy flux has components 
T n and which satisfy 20 

[ fid dJ 7 ^ 

oj a J dm a I a + V k - T a = u a j dm a I a + -—{nfZ) + -^=0. (AlO) 

Care must be taken in evaluating eq.A8. For finite values of n a and taking the limit (p, q) — > (po,qo) 
gives J dm a I a — 0. However, in the vicinity of n a = or fi a = eq.A8 takes on the indeterminate form 
0/0, since X — > in the numerator and either n a or [i a goes to zero in the denominator. Following Zakharov 
et.al., we evaluate this expression using the delta function representation lim £ ^ e|a ; | e_1 = 28(x) to find 

dm a I a = -4N$ I — -z-(po,qo)8(n a ) H -s~(Po, Qo)S(^ a ) ) ■ (All) 

VMa dp n a dq J 

Hence, the flux can only have a source for n = or \i = 0. Plugging eq.All into eq.AlO, using ui = 2/i, and 
integrating over n or /x, respectively, gives the two components of the flux 

dJ 

T n = 8JV?n-V _1 ^-(Po,«b) 

?» = 8N^n- 2 ^-(p ,q a ). (A12) 

We now estimate the dimensionless flux integrals, verifying that they converge and give have the correct 
sign. We do this by breaking the integration up into two regimes: large daughter mode wavenumber (xp >• 1) 



20 The assumption of local energy transfer is later shown to be valid by verifying that the flux integrals converge. 
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and small daughter mode wavenumber {xp <C 1). Although our expansions for the integrand are technically 
only valid in the respective limits, we extend them all the way to 1,3 ~ 1. 

When the derivatives with respect to p and q arc taken in eq.A9, only the last parenthesis need be 
differentiated since it gives zero for (p, q) = (po,qo)- The coupling coefficients for the two limits are given in 
appendix B. The resulting expressions are 

&X f°° 128 

0^(Po> 90 ) - 2 3/2 J dx fj Xp 2 lnx fj + — J dxpx 1 / 2 (lnx fj + 1) 

= 2 3 / 2 + S-3.7.« n (A13) 



and 
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-7^(Po,qo) = -2 3/2 ln2y dx p xj 2 + —3- J dxpx^lnxp 

= -2 3 / 2 ln2- — = -3.8=-a„. (A14) 

We note that the contribution to T 11 ^^) from both large and small xp arc positive (negative). 

The energy flux is toward larger n and smaller p,. Wc now restore the sign of p in order to have an 
expression valid for either sign. The final result for the energy flux is then 

T» = -8a^N 2 n- 2 -^. (A15) 
iMl 

We also note the final answer for the occupation number 

N = N n- 4 fi- 1/2 (A16) 
where N a is related to the energy flux by eq.A15. 



B. analytic estimate of the maximum coupling coefficient 

Here we give a detailed analytic calculation of the maximum coupling coefficient. We make the following 
approximations: (1) n 3> m, the WKB limit; (2) k§ <C &r, which follows from (1); (3) short wavelength 
daughter modes with pp ~ p 1 ~ — fj, a /2; (4) an n = 1 polytropic background star. For convenience, we 
use units with r* = M = 1, and we use the normalization condition £Wit = 0.5Mr^O 2 . In these units, the 
density profile near the surface is p = poz, where the central density is po = 7r/4. 

We decompose tp in eq.42 as a sum of plane waves as 



, / \-l/2 

= T^rW ( — sin dl sin 62 ) ex P( i Xs 1 s 2 ) (Bl) 



where the WKB phase is 



= p(s 1 e 1 +s 2 6 2 ) +m^ + a(si + s 2 ). (B2) 



In the small p limit, we can write 9 2 — tt/2 — \p\e where — 1 < e < 1. Hence the effective wavenumber in the 
6*2 direction is p\p,\ = ir\k\. From eq.29 we find the displacement vector 
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where 

e R = (k R + iqk^) ~ -ip?k R 

1 — q z 
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£</> = 1 _ g2 ( k 4> ~ ~ > £fl 

e z = ifc z ~ e . (B4) 
The wavenumber is defined by fc,s lS2 = dx Sl s 2 /dx with components found from eq.34 to be 

k R = P{1 ~ ( SlXl{1 _ ^2)1/2 _ S2X2{1 _ X \f/A 

m to(1 — p 2 ) 1 ^ 2 



fc 2 = -5i52 I" 1/2 fc fl . (B5) 

The factor (x 2 — x 2 )/|^i| is the volume element for the bi-spheroidal coordinates, and becomes zero at 
coordinate singularities. In the limit z«l and fi <C 1, it has the simple form 

(* 2 -* 2 ) 2 * ( y 2 - ?y 2 ) 2 + A 2 (B6) 

where y = cos#, yo — (p 2 — 25) 1 / 2 , and A = (Sp^z) 1 ^ 2 . This formula shows that there is a narrow peak for 
the wavenumber (for z < p 2 /2) near the singular point (r, cos#) = (l,±|/z|); there is no pronounced peak 
for z > /U 2 /2 and the integrand decreases strongly. 

Plugging the WKB travelling wave forms into eq.49 and keeping only the largest terms in the WKB limit 
gives 
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where the matrix element is defined by 

M a p 7 — yt 2 k a ■ epk a ■ e 7 + p 2 kp • e 7 kp ■ e a + p 2 1 k 1 ■ e a fc 7 • ep (B8) 

and p a is the cutoff version of the WKB envelope. 

To evaluate this expression, we first note that the integral is rapidly oscillating unless the conservation 
rules of eq. 50 are satisfied. These relations are written 

SalPa + Spipp + S 7 ip 7 = 

s a 2\k a \ + sp 2 \kp\ + s 72 |/c 7 | = 

m a + rap + to 7 = 0. (B9) 

In the limit of large daughter mode wavenumber, we find pp ~ p 7 , s 7 i = — sp\, |fc 7 | ~ \kp\, s 72 = — sp2, and 
to 7 ~ — TO/3; m other words, the wavevectors of the daughter modes are equal in magnitude and opposite in 
direction. In this limit, the remaining spatially constant phase factor is 

exp[i(x„ + xp + Xj)} = i-^ s ^+ s ^+^- s ^ s ^+ s ^ (BIO) 
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where 5 = for an even parity mode and <5 = 1 for an odd parity mode. Using the incompressibility condition 
and momentum conservation leads to the simplification kp ■ e 7 = — fc Q • e 7 ~ k a -ep. Plugging these relations 
in gives 

M afil = ^p 2 a k a ■ ep {-2k a -£p + kp- e a ) . (Bll) 
In the limit of small p and k$ <C k R we find 

k a -ep ~ -i-^k aR kp R + ik az kp z 

kp-e a ~ -ip 2 a k aR kp R + ik az kp z . (B12) 
Performing the spin sums, we find zero for net odd parity, while for even parity we get 



6 2 2 

V M a ,-[«a(»al+»a2) + («/9-^)(^l+^2)] _ MqPqP/? 

Jw a/?7' ^ / 2 _ „2 \2rfV 2 _ T 2 \2 



x - a£ 2 ) + X L(1 - x 'i)} - ^2) + ^2(! - 4i)} 

+ 4s a ia; a 2a; / jia: / j2{(l-a:^ 1 )(l-a^ 2 )(l-4i)(l-42)} 1/2 )- ( B13 ) 

So far, we have 

4£l 2 Po V'OaV'gfl 6 2 2 Pa J 

= 2£ unit (2*)« W(iJ = 4^ (B14) 



where the dimensionless integral is 

/2 

-1/2 



, 1/2 



J = [ d 3 x— (sin 0^ sin q2 sin 2 0«i sin 2 0,32)" 

7 PO \ PaPg / 



2 I ^.ui-.ai- „ pi „.„ „ - 2 _ T 2 \2 

*F/3 / ^al x a2,l l x /31 x p2> 

x ({^(i - x 2 2 ) + z 2 2 (i - x 2 !)} {4i(i - 4 2 ) + 4 2 (i - 4i)} 

+ 4x Ql x Q 2 x /3iX /32 {(l-2;L)( 1 -a ; L)( 1 -4i)( 1 "42)} 1/2 ) • ( B15 ) 

To evaluate the dimensionless integral we first note that if one attempted to take the p — > limit, as 
for the mode energy, one would find an integral 

17 a /„ (B16) 
which implies a linear divergence at R = 1. The divergence implies a strong dependence of the integral on 
some characteristic lengthscale near the surface. 

To proceed with a more detailed calculation, first notice that the integrand is sharply peaked near the 
surface. Below the daughter mode turning point zp, a factor of density cancels out. The parent mode 
WKB amplitude is largest for above the turning point z < z a , where the cutoff density for the parent mode 
Pa/ 'po ~ z a . Near the daughter mode singular point cos 9 = y = \pp\, the bi-spheroidal coordinates become 

X a i — \Pa\i x a2 

— IMjsU xpi ^ x fj2 ^\pp\- A similar expression can be found near the parent mode singular 
point. The integral J then becomes 

dy 



= Wnz- 1 / 2 ["" dz-^ + 4nz^ 2 ["" dz-^ 
4V2tt 2 f Ba z 4V2tt 2 r<* z 16V2vr 2 



(j/2- y 2 o)2+A 2 



4V2tt 2 r* a J_ W2ir 2 p _z_ _ 

~~ -1/2 2 /- ?l/2 + -1/2 , / *T/2 ~ 



Hi 



(B17) 
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The final result for the coupling coefficient for nearly identical, short lengthscale daughter modes is then 



— — 



p a 16\/2tt 2 

4^V2 M 2 



1.02m^ 2 . 



(B18) 



In appendix A when finding the cascade solution, one also needs the coupling coefficient in the limit 
that one of the daughter modes has a small wavenumbcr. In section 4.1 we found that the approximate 
solutions in the limit n 7 <C n a ~ n 7 are: n 7 ~ n a + np, fc 7 ~ k a , hp ~ k a (rii3/n a ) 2 , /i 7 ~ — /i Q , and 
fij3 ~ H a (np/n a ) <C Since the method is exactly the same as the detailed calculation already given, we 
merely quote the answer for the maximum coupling coefficient to be 



Here So — 1 if modes a and 7 have the opposite parity, and So — otherwise. Similarly, 5s = 1 if modes a 
and 7 have the same parity, and 5s = otherwise. 




(B19) 
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